© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration 2013APlanck Collaboration I 2014), describes new measurements of the cosmic infrared background (CIB) anisotropy power spectrum and bispectrum, and their use in constraining the cosmic evolution of the star formation density and the luminous-dark matter bias.

The relic emission from galaxies formed throughout cosmic history appears as a diffuse, cosmological background. The CIB is the far-infrared part of this emission and it contains about half of its total energy (Dole et al. 2006). Produced by the stellar-heated dust within galaxies, the CIB carries a wealth of information about the process of star formation. Because dusty, star-forming galaxies at high redshift are extremely difficult to detect individually (e.g., Blain et al. 1998; Lagache et al. 2003; Dole et al. 2004; Fernandez-Conde et al. 2008; Nguyen et al. 2010), the CIB represents an exceptional tool for studying these objects and for tracing their overall distribution (Knox et al. 2001). The anisotropies detected in this background light trace the large-scale distribution of star-forming galaxies and, to some extent, the underlying distribution of the dark matter haloes in which galaxies reside. The CIB is thus a direct probe of the interplay between baryons and dark matter throughout cosmic time.

The CIB has a redshift depth that complements current optical or near infrared measurements. This characteristic can be used to explore the early build-up and evolution of galaxies, one of the main frontiers in cosmology. Indeed the hope is to be able to use CIB anisotropies to improve our understanding of early gas accretion and star formation, and to assess the impact of galaxies on reionization. As dusty star-forming galaxies start to be found up to very high redshift (e.g., z = 6.34, Riechers et al. 2013), this objective may be reachable, even if quantifying the z ≲ 6 contribution to CIB anisotropy measurements to isolate the high-redshift part will be very challenging. However, as a start, CIB anisotropies can be used to measure the cosmic evolution of the star formation rate density (SFRD) up to z ≃ 6. Quantifying the SFRD at high redshift (z> 2.5) is a challenging endeavour. Currently, most of the measurements rely on the UV light emerging from the high-redshift galaxies themselves (e.g., Bouwens et al. 2009; Cucciati et al. 2012). To estimate their contribution to the total SFRD one needs to apply the proper conversion between the observed UV rest-frame luminosity and the ongoing SFR. This conversion factor depends on the physical properties of the stellar population (initial mass function, metallicities and ages) and on the amount of dust extinction, and is thus rather uncertain. Despite the significant amount of effort aimed at better understanding the UV-continuum slope distribution at high redshift, this remains one of the main limitation to SFRD measurements. The uncertainty on this conversion sometimes leads to significant revision of the SFRD (e.g., Behroozi et al. 2013b; Bouwens et al. 2012; Castellano et al. 2012). Remarkably, the estimates are now routinely made up to z ~ 8 (e.g., Oesch et al. 2012a), and have even been pushed up to z ~ 10 (Oesch et al. 2012b). One of the key questions, is how to quantify the contribution of the dusty, star-forming galaxies to the SFRD at high redshift. Since it is mostly impossible to account for this contribution on the basis of optical/near-IR surveys, the best approach is to use the dusty galaxy luminosity function measurements. However, such measurements at high redshift are challenging with the current data, and this is where the CIB anisotropies, with their unmatched redshift depth, come into play. The SFRD from dusty, star-forming galaxies can be determined from their mean emissivity per comoving unit volume, as derived from CIB anisotropy modelling.

The way galaxies populate dark matter haloes is another ingredient that enters into the CIB anisotropy modelling. In particular, the galaxy bias – the relationship between the spatial distribution of galaxies and the underlying dark matter density field – is a result of the varied physics of galaxy formation which can cause the spatial distribution of visible baryons to differ from that of dark matter. If galaxy formation is mainly determined by local physical processes (such as hydrodynamics), the galaxy bias is then approximately constant on large scales (Coles 1993), and the galaxy density fluctuations are thus proportional to those of the dark matter. The proportionality coefficient here is usually called as the linear bias factor, b. Its dependence on the luminosity, morphology, mass, and redshift of galaxies provides important clues to how galaxies are formed. However, the linear biasing parameter is at best a crude approximation, since the true bias is likely to be nontrivial, i.e., non-linear and scale dependent, especially at high redshift. At high redshift, the biasing becomes more pronounced, as predicted by theory (e.g., Kaiser 1986; Mo & White 1996; Wechsler et al. 1998), and confirmed by the strong clustering of dusty star-forming galaxies (e.g., Steidel et al. 1998; Blain et al. 2004; Cooray et al. 2010). CIB anisotropies can be used to constrain the biasing scheme for dusty star-forming galaxies, which is crucial for understanding the process and history of galaxy formation.

However, measuring the CIB anisotropies is not easy. First, the instrument systematics, pipeline transfer function and beams have to be very well understood and measured. One can take advantage of recent experiments such as Hershel and Planck, for which diffuse emission is measured with better accuracy than their IRAS and Spitzer predecessors. Second, extracting the CIB requires a very accurate component separation. Galactic dust, CMB anisotropies, emission from galaxy clusters through the thermal Sunyaev Zeldovich (tSZ) effect, and point sources all have a part to play. In clean regions of the sky, Galactic dust dominates for multipoles ≲200. This has a steep power spectrum (with a slope of about − 2.8) and exhibits spatial temperature variations, and thus spectral energy distribution (SED) spatial variations. Distinguishing Galactic from extragalactic dust is very difficult, as their SEDs are quite similar, and both their spatial and spectral variations do not exhibit any particular features. Currently, the best approach is to rely on a Galactic template; taking another frequency is not recommended as CIB anisotropies also contribute. Taking a gas tracer as a spatial template is the best one can do, even if it has the drawback of not tracing the dust in all interstellar medium phases. For this purpose the CMB is very problematic at lower HFI frequencies, since its power spectrum is about 5000 and 500 times higher at ℓ = 100 than the CIB at 143 and 217 GHz, respectively. For Planck  at 217 GHz  the CMB dominates the CIB for all ℓ; at 353 GHz  it dominates for ℓ < 1000; and at 545 GHz  its power is 25 times lower at ℓ = 100 than the CIB. Any CMB template (taken from low-frequency data or from complex component separation algorithms) will be contaminated by residual foreground emission that will have to be corrected for. At Planck frequencies ν> 200  GHz  the tSZ effect can be safely ignored, being close to zero at 217  GHz and 100 times below the CIB power spectrum at 353  GHz. However tSZ contamination can come from the use of a CMB template that contains residual tSZ power. At 100  GHz, based on the tSZ power spectrum measured in Planck Collaboration 2013UPlanck Collaboration XXI (2014) and the CIB model developed in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), we estimate the tSZ power spectrum to be 10 times higher than the CIB. The correction of any tSZ contamination will be made difficult by the intrinsic correlation between the CIB and tSZ signals (Addison et al. 2012; Reichardt et al. 2012). Finally, as bright point sources will put extra power at all scales in the power spectrum, point sources need to be carefully masked up to a well-controlled flux density level. This step is complicated by the extragalactic source confusion that limits the depth of source detection for current far-infrared and submillimetre space missions.

The pioneering studies in CIB anisotropy measurement with Herschel-SPIRE (Amblard et al. 2011) and Planck-HFI (Planck Collaboration 2011RPlanck Collaboration XVIII 2011) are now extended in Viero et al. (2013b) for the former, and in this paper for the latter. The new Planck measurements benefit from larger areas, lower instrument systematics, and better component separation. They are not limited to auto-power spectra but also include frequency cross-spectra, from 143 to 857  GHz (and 3000 GHz with IRAS), and extend to bispectra at 217, 353, and 857  GHz. These more accurate measurements pose new challenges for the modelling of CIB anisotropies.

Our paper is organized as follows. We present in Sect. 2 the data we are using and the field selection. Section 3 is dedicated to the removal of the background CMB and foreground Galactic dust. We detail in Sect. 4 how we estimate the power spectrum and bispectrum of the residual maps, and their bias and errors. In the same section, results on CIB power spectra and bispectra are presented. In Sect. 5, we describe our modelling and show the constraints obtained on the SFRD, and the clustering of high-redshift, dusty galaxies. In Sect. 6, we discuss the 143  GHz anisotropies, the frequency decoherence, the comparison of our measurements with previous determinations, the SFRD constraints, and the CIB non-Gaussianity. We conclude in Sect. 7. The appendices give some details about the Hi data used to remove Galactic dust (Appendix A), the CIB anisotropy modelling (Appendices B and C), and also present the power spectra and bispectra tables (Appendix D).

Throughout the paper, we adopt the standard ΛCDM cosmological model as our fiducial background cosmology, with parameter values derived from the best-fit model of the CMB power spectrum measured by Planck (Planck Collaboration 2013PPlanck Collaboration XVI 2014): { Ω_m,Ω_Λ,Ω_bh²,σ₈,h,n_s } = { 0.3175,0.6825,0.022068,0.8344,0.6711,0.9624 }. We also adopt a Salpeter initial mass function (IMF).

2. Data sets and fields

2.1. Planck HFI data

We used Planck channel maps from the six HFI frequencies: 100, 143, 217, 353, 545, and 857  GHz. These are N_side = 2048HEALPix² maps (Górski et al. 2005), corresponding to a pixel size of 1.72′. We made use of the first public release of HFI data that corresponds to temperature observations for the nominal Planck mission. The characteristics of the maps and how they were created are described in detail in the two HFI data processing and calibration papers (Planck Collaboration 2013FPlanck Collaboration VI 2014; Planck Collaboration 2013HPlanck Collaboration VIII 2014). At 857, 545, and 353 GHz, we use the zodiacal light subtracted maps (Planck Collaboration 2013NPlanck Collaboration XIV 2014). Some relevant numbers for the CIB analysis are given in Table 1.

Maps are given in units either of MJy sr^-1 (with the photometric convention νI_ν = const.)³ or K_CMB, the conversion between the two can be exactly computed knowing the bandpass filters. The mean coefficients used to convert frequency maps in K_CMB units to MJy sr^-1 are computed using noise-weighted band-average spectral transmissions. The mean conversion factors are given in Table 1. The map-making routines do not average individual detector maps, but instead combine individual detector data, weighted by the noise estimate, to produce single-frequency channel maps. As portions of the sky are integrated for different times by different detectors, the relative contribution of a given detector to a channel-average map varies for different map pixels. The effects of this change on the channel-average transmission spectra is very small, being of the order of 0.05% for the nominal survey coverage (Planck Collaboration 2013IPlanck Collaboration IX 2014; Planck Collaboration 2013).

As in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), the instrument noise power spectrum on the small extragalactic fields (see Sect. 2.3) is estimated with the jack-knife difference maps, which are built using the first and second halves of each pointing period (a half-pointing period is of the order of 20 min). The half-ring maps give an estimate of the noise that is biased low (by a couple of percent) due to small correlations induced by the way the timelines have been deglitched (Planck Collaboration 2013JPlanck Collaboration X 2014). However, as discussed in Planck Collaboration 2013FPlanck Collaboration VI (2014), the bias is significant only at very high multipoles. As we stop our CIB power spectra measurements at ℓ = 3000, we can safely ignore this bias. For the larger GASS field (see Table 2), we directly compute the cross-spectra between the two half-maps to get rid of the noise (assuming that the noise is uncorrelated between the two half-maps).

The effective beam window functions, b_ℓ, are determined from planet observations, folding in the Planck scanning strategy, as described in Planck Collaboration 2013GPlanck Collaboration VII (2014). Because of the non-circular beam shape, the detector combinations, and the Planck scanning strategy, the effective b_ℓ(ν) of each channel map applicable to the smallest patches considered here varies across the sky. These variations are, however, less than 1% at ℓ = 2000 and average out when considering many different patches, or larger sky areas. We will therefore ignore them and consider a single b_ℓ(ν) for each frequency channel. Because of their experimental determination, the b_ℓ(ν) are still affected by systematic uncertainties, which can be represented by a small set of orthogonal eigen-modes and whose relative standard deviation can reach 0.5% at ℓ = 2000 for ν = 857 GHz. This uncertainty is accounted for in the error budget.

2.2. IRIS data

Our analysis uses far-infrared data at 3000 GHz (100μm) from IRAS (IRIS, Miville-Deschênes & Lagache 2005). During its 10-month operation period, IRAS made two surveys of 98% of the sky and a third one of 75% of the sky. Each survey, called an HCON for Hours CONfirmation, consisted of two coverages of the sky separated by up to 36 h. The first two HCONs (HCON-1 and HCON-2) were carried out concurrently, while the third survey (HCON-3) began after the first two were completed. Due to exhaustion of the liquid helium supply, the third HCON could not be completed. We use the HCON difference maps to estimate the instrument noise power spectrum (just as we use the half-pointing-period maps to estimate the instrument noise power spectrum for HFI). As for HFI at high frequencies, the IRIS 3000 GHz map is given in MJy sr${-1}^{}$ with the photometric convention νI_ν = const.

The beam at 3000 GHz is not as well characterized as the HFI beams. The IRIS effective FWHM is about 4.3′ (Miville-Deschênes & Lagache 2005). We estimate the FWHM uncertainty using different measurements of the point spread function (PSF) coming from selected point sources used to study the PSF in the “ISSA explanatory supplement”, and the power spectrum analysis from Miville-Deschênes et al. (2002). The dispersion between those estimates gives an uncertainty of 0.5′ on the FWHM.

2.3. Extragalactic fields with high angular resolution HI data

Following the successful approach of Planck Collaboration 2011RPlanck Collaboration XVIII (2011), we do not remove Galactic dust by fitting for a power-law power spectrum at large angular scales, but rather use an independent, external tracer of diffuse dust emission, the Hi gas. From 100μm to 1 mm, 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; Planck Collaboration 2011XPlanck Collaboration XXIV 2011). 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.

Although Planck is an all-sky survey, we restricted our CIB anisotropy measurements to a few fields at high Galactic latitude, where Hi data at an angular resolution close to that of HFI are available. The 21-cm Hi spectra used here were obtained with: (1) the Parkes 64-m telescope; (2) the Effelsberg 100-m radio telescope; and (3) the 100-m Green Bank telescope (GBT). Field characteristics are given in Table 2. Further details on the Hi data reduction and field selection are given in Appendix A. The HEALPix HFI maps were reprojected onto the small Hi GBT and EBHIS maps by binning the original HEALPix data into Hi map pixels (Sanson-Flamsteed, or “SFL” projection with pixel size of 3.5′ for all fields). An average of slightly more than four HEALPix pixels were averaged for each small map pixel. For the GASS field, the HFI data were convolved to the Hi angular resolution (16.2′), and then degraded to N_side = 512.

We have 11 fields (one EBHIS, nine GBT and one GASS). The total area used to compute the CIB power spectrum is about 2240 deg², 16 times larger than in Planck Collaboration 2011RPlanck Collaboration XVIII (2011).

2.4. Point sources: flux cut and masks

We use the Planck Catalogue of Compact Sources (PCCS, Planck Collaboration 2013ZCPlanck Collaboration XXIX 2014) to identify point sources with signal-to-noise ratio (S/N) greater or equal to 5 in the maps, and we create point-source masks at each frequency. We mask out a circular area of 3σ radius around each source (where σ = FWHM/ 2.35). The point sources to be removed have flux densities above a chosen threshold. The threshold is determined using the number counts dN/dS on the cleanest 30% of the sky and by measuring the flux density at which we observe a departure from a Euclidean power law; this departure is a proxy for measuring the flux density regime where the incompleteness starts to be measurable. In practice, this departure often corresponds to about 80% of completeness (e.g., Planck Collaboration Int GPlanck Collaboration Int. VII 2013). The uncertainty on flux densities comes from the form of the dN/dS slope. Flux density cuts and all-sky effective beam widths are given in Table 1. For IRIS at 3000 GHz, we use the IRAS faint source catalogue (Moshir 1992) and we mask all sources with S₃₀₀₀> 1 Jy.

3. Extracting CIB from Planck-HFI and IRIS maps

One of the most difficult steps in extracting the CIB is the removal of the Galactic dust and the CMB. CMB anisotropies contribute significantly to the total HFI map variance, in all channels at frequencies up to and including 353  GHz. Galactic dust contributes at all frequencies, and is dominant at the higher frequencies. One approach is to keep all the components and search for the best-fit model of the CIB in a likelihood approach, accounting for both CMB and dust (e.g., as a power law). This is the philosophy of the method developed to extract the cosmological parameters from the lower frequency data analysis (Planck Collaboration 2013OPlanck Collaboration XV 2014; Planck Collaboration 2013PPlanck Collaboration XVI 2014). However, with such an approach, the complexity of the likelihood, together with the number of parameters and their degeneracies prevent the use of advanced models for the clustered CIB (models beyond a simple power law). We thus decided to use another approach, based on template removal. To remove the CMB in the fields retained for our analysis, we used a simple subtraction technique (as in Planck Collaboration 2011RPlanck Collaboration XVIII 2011). This method enables us to reliably evaluate CMB and foreground component residuals, as well as noise contamination, and to easily propagate errors (for example cross-calibration errors). It also guarantees that high-frequency CIB anisotropy signals will not leak into lower frequency, CMB-free maps. For Galactic dust, the present work focusses on very clean regions of the sky, for which Galactic foregrounds can be safely monitored using ancillary Hi observations.

We describe in this section the removal of the two components. Some corrections are made at the map level, while others can only be done at the power spectrum level.

3.1. CMB removal

3.1.1. A low-frequency map as a CMB template

The extraction of CIB anisotropies at low frequency is strongly limited by our ability to separate the CIB from the CMB. As a matter of fact, at multipole ℓ ~ 100, the CIB anisotropy power spectrum represents about 0.2% and 0.04% of the CMB power spectrum at 217 and 143  GHz, respectively. We decided in this paper to use the HFI lowest frequency channel (100  GHz) as a CMB template. This has the advantage of being an “internal” template, meaning its noise, data reduction processing, photometric calibration, and beam are all well known. It also has an angular resolution close to that of the higher frequency HFI channels. Following Planck Collaboration 2011RPlanck Collaboration XVIII (2011), we removed the CMB contamination in the maps at ν ≤ 353 GHz. At 545  GHz, because the CMB power represents less than 5% of the CIB power, we removed the CMB in harmonic space.

Using the lowest frequency channel as a CMB template ensures the lowest CIB contamination, since the CIB SED is decreasing as ν^3.4 (Gispert et al. 2000). Contrary to what was done in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), where the 143  GHz channel was used as a CMB template, we used the HFI 100  GHz map here. This is a good compromise between being a low frequency template, and having an angular resolution close to the higher frequency HFI channels. Note that using the 100  GHz channel was not feasible in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), because the data were too noisy, and the field area too small (140 versus 2240 deg²).

As detailed in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), we applied a Wiener filter to the 100  GHz map, designed to minimize the contamination of the CMB template by instrument noise. Errors in relative photometric calibration (between channels) are accounted for in the processing, as detailed in Sects. 4.1 and 4.2.

Following Planck Collaboration 2011RPlanck Collaboration XVIII (2011), two corrections have to be applied to the measured CIB power spectra when using such a CMB template. First we need to remove the extra instrument noise that has been introduced by the CMB removal. This is done through: ${\mathit{N}}_{\mathit{\ell }}^{\mathrm{CMBres}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{100}}\mathrm{\right)}\mathrm{\times }{\mathit{w}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\times }{\left(\frac{{\mathit{b}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}}{{\mathit{b}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{100}}\mathrm{\right)}}\right)}^{\mathrm{2}}\mathit{,}$(1)with ν equal to 143, 217 or 353 GHz. b_ℓ(ν) is the beam window function, w_ℓ is the Wiener filter, and N_ℓ(ν₁₀₀) is the noise power spectrum of the 100  GHz map. It is computed in the same way as in the other frequency channels, using the half-pointing period maps. Second, owing to the lower angular resolution of the 100  GHz channel compared to 143, and 217 and 353  GHz, we also have to remove the CMB contribution that is left close to the angular resolution of the 143, 217 and 353 GHz channels: ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMBres}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{\times }{\mathit{F}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\times }{\mathit{b}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{\times }{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{w}}_{\mathit{\ell }}\mathrm{)}}^{\mathrm{2}}\mathit{,}$(2)with F_ℓ being the pixel and reprojection transfer function (detailed in Sect. 4.1).

Note that Eqs. (1) and (2) are the corrections for the auto-power spectra. They are easily transposable to cross-power spectra.

3.1.2. CMB template contamination to the CIB

The low-frequency channel CMB template has the disadvantage of being contaminated by the CIB and the tSZ effect (the Galactic dust is removed using our dust model detailed in the following section, and IR and radio point sources are masked). Indeed the cross-spectra between the estimated CIB maps at ν and ν′ involve the a_ℓm product $\begin{array}{ccc}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{\nu }}\mathrm{\times }{\mathit{a}}_{\mathit{\ell m}}^{{\mathit{\nu }}^{\mathrm{\prime }}\mathrm{\ast }}& \mathrm{=}& \left\{{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{CIB}\mathit{,\nu }}\mathrm{+}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{SZ}\mathit{,\nu }}\mathrm{-}{\mathit{w}}_{\mathit{\nu }}\mathrm{(}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{CIB}\mathit{,}\mathrm{100}}\mathrm{+}{{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{SZ}\mathit{,}\mathrm{100}}}^{\mathrm{)}}\right\}\\ & & \end{array}$(3)where w_ν is either zero if ν ≥ 545, or the Wiener filter applied to the 100 GHz map (w_ℓ) if ν ≤ 353.

Besides the signal ${\mathit{C}}_{\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ that we want to measure, Eq. (3) involves three additional contributions that we have to correct for: tSZ × tSZ; spurious CIB × CIB; and CIB × tSZ correlations. We discuss each of them in this section.

CIB × CIB spurious correlations –

From Eq. (3), the CIB × CIB spurious contribution reads ${\mathit{C}}_{\mathrm{CIBcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}\mathrm{-}{\mathit{w}}_{\mathit{\nu }}{\mathit{C}}_{\mathrm{CIB}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{-}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{C}}_{\mathrm{CIB}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{\nu }}\mathrm{+}{\mathit{w}}_{\mathit{\nu }}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{C}}_{\mathrm{CIB}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}}\mathit{.}$(4)Using a model of CIB anisotropies, we can compute ${\mathit{C}}_{\mathrm{CIBcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$. This correction includes both the shot noise and the clustered CIB anisotropies. The correction will be taken into account when searching for the best-fit CIB model: for each realization of our CIB anisotropy models (detailed in Sects. 5.4 and 5.5) we will compute ${\mathit{C}}_{\mathrm{CIBcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$.

To illustrate the order of magnitude of the correction, we can use the model constructed in Planck Collaboration 2011RPlanck Collaboration XVIII (2011) to fit the early Planck CIB measurements, and the shot-noise levels given in Sect. 5.1, and compute ${\mathit{C}}_{\mathrm{CIBcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ when subtracting the 143 GHz map instead of the 100 GHz map as a CMB template. Note that this comparison is only indicative, since this CIB model at 143 and 100  GHz is an extrapolation of the 217 × 217 power spectrum, using the clustering parameters of the 217  GHz best-fit model, and the emissivities computed using the same empirical model of galaxy evolution. We show this comparison for the 217 × 217 power spectrum in Fig. 1. We see from this figure that the CIB obtained using the 143  GHz map as a CMB template is about 1.6 times lower than the CIB obtained using the 100  GHz map. Applying the correction (Eq. (4)) largely decreases this discrepancy (compare the two dashed lines⁵). Note that the correction is also non-negligible when using the 100  GHz map: it is a factor 1.15 for 50 <ℓ < 700. This justifies applying the correction systematically when fitting for our best CIB models.

We can extend the check on the impact of the choice of CMB map using the cross-correlation between the CIB and the distribution of dark matter, via the lensing effect on the CMB (Planck Collaboration 2013RPlanck Collaboration XVIII 2014). We specifically compared the cross-correlation for the CIB at 217  GHz obtained using our two Wiener-filtered 100 and 143  GHz maps. We reached the same conclusion as before, and retrieved the same underestimate of the CIB when using the Wiener-filtered 143  GHz rather than 100  GHz CMB map.

Finally, we note that an alternative method of removing the CMB contamination was extensively tested; this was based on an internal linear combination of frequency maps, combined (or not) with a needlet analysis. However, such CMB maps are often not suited to our purposes because they are, among other problems, contaminated by the CIB that has leaked from the high-frequency channels that are used in the component separation process. We compare the CIB obtained using the 100  GHz and SMICA CMB maps (Planck Collaboration 2013LPlanck Collaboration XII 2014) in Fig. 1. We see that they are compatible within 1σ (red-dashed line and orange points). Again this comparison is only indicative, since the SMICA CMB map main contaminants are SZ, shot noise from point sources, and CIB⁶; those contaminants have not been corrected for here, because the correction would need fully realistic simulations.

tSZ×tSZ correlation

From Eq. (3), the tSZ×tSZ spurious contribution reads ${\mathit{C}}_{\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}{\mathit{C}}_{\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{+}{\mathit{w}}_{\mathit{\nu }}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{C}}_{\mathrm{SZ}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}}\mathrm{-}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{C}}_{\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}}\mathrm{-}{\mathit{w}}_{\mathit{\nu }}{\mathit{C}}_{\mathrm{SZ}}^{{\mathit{\nu }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}}\mathit{.}$(5)We can make explicit the frequency dependence of the tSZ and write $\mathrm{\Delta }{\mathit{T}}_{\mathrm{SZ}}^{\mathit{\nu }}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}{\mathit{g}}_{\mathit{\nu }}\sum _{\mathit{\ell m}}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{SZ}}\mathit{,}$(6)where g_ν is the conversion factor from the tSZ Compton parameter y to CMB temperature units. Hence ${\mathit{C}}_{\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{\nu }}\mathrm{=}{\mathit{g}}_{\mathit{\nu }}^{\mathrm{2}}{\mathit{C}}_{\mathrm{SZ}}\mathit{,}$(7)Eq. (5) then becomes ${\mathit{C}}_{\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}{\mathit{C}}_{\mathrm{SZ}}\mathrm{\{}{\mathit{g}}_{\mathit{\nu }}{\mathit{g}}_{{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{+}{\mathit{w}}_{\mathit{\nu }}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{g}}_{\mathrm{100}}^{\mathrm{2}}\mathrm{-}{\mathit{g}}_{\mathrm{100}}\mathrm{(}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{g}}_{\mathit{\nu }}\mathrm{+}{\mathit{w}}_{\mathit{\nu }}{\mathit{g}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathrm{\left)}}^{\mathrm{\right\}}}\mathit{.}$(8)We compute C_SZcorr using the tSZ power spectrum, C_SZ, and the conversion factors, g_ν, given in Planck Collaboration 2013UPlanck Collaboration XXI (2014). The uncertainty on C_SZ is about 10%.

tSZ×CIB correlation

From Eq. (3), the tSZ×CIB spurious contribution reads $\begin{array}{ccc}{\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}& \mathrm{=}& {\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{+}{\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{{\mathit{\nu }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{\nu }}\mathrm{-}{\mathit{w}}_{\mathit{\nu }}{\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{-}{\mathit{w}}_{{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{\nu }}\\ & & \end{array}$(9)where ${\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ is the notation for the cross-spectrum between CIB(ν) and SZ(ν′). This correction is highly dependent on the model used to compute the cross-correlation between tSZ and CIB; we use the model from Addison et al. (2012). We make the assumption that ${\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}\mathrm{⟨}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{CIB}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{SZ}\mathrm{\ast }}\mathrm{\left(}{\mathit{\nu }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{g}}_{{\mathit{\nu }}^{\mathrm{\prime }}}\mathit{\phi }\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZ}}\mathit{,}$(10)where φ(ν) is the amplitude of the power spectrum of the CIB correlated with the tSZ, C_{CIB × SZ}, taken from Addison et al. (2012). This paper also provides examples of cross-spectra, with a reference frequency at 150 GHz. We use this reference frequency and power spectrum ratios to compute ${\mathit{C}}_{\mathrm{SZ}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ following ${\mathit{C}}_{\mathrm{SZ}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}\frac{\mathit{\phi }\mathrm{\left(}{\mathit{\nu }}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\phi }\mathrm{\left(}\mathrm{150}\mathrm{\right)}}\frac{\mathit{g}\mathit{\nu }}{{\mathit{g}}_{\mathrm{150}}}{\mathit{C}}_{\mathrm{SZ}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{CIB}}^{\mathrm{150}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{150}}\mathit{.}$(11)We show in Fig. 2 the measured CIB, ${\mathit{C}}_{\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathrm{measured}\mathrm{\right)}$, together with the corrected one: ${\mathit{C}}_{\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}{\mathit{C}}_{\mathrm{CIB}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathrm{measured}\mathrm{\right)}\mathrm{-}\mathrm{\{}{\mathrm{C}}_{\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{+}{{\mathrm{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}}^{\mathrm{\}}}\mathit{.}$(12)We also show the ratio of the corrections to the measured CIB, and give the values of the corrections at 217 and 353 GHz in Table 3. We see that tSZ contamination is the highest for the 217 × 217 combination, with a contamination of the order of 15%. It is less than 5% and 1% for 217 × 353 and 353 × 353, respectively. The tSZ power spectrum measured by Planck Collaboration 2013UPlanck Collaboration XXI (2014) is uncertain by 10%, while tSZ×CIB is uncertain by a factor of two (Addison et al. 2012). Hopefully, where the contamination is important (i.e., 217 × 217), the dominant contribution to the CIB comes from ${\mathit{C}}_{\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$. Hereafter, we therefore apply the correction coming from the tSZ contamination to the measured CIB, and add the uncertainties of the correction quadratically to the CIB uncertainties.

When cross-correlating maps at 353  GHz and above 545  GHz, the correction linked to the tSZ contamination is dominated by the term ${\mathit{C}}_{\mathrm{CIB}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$, which is highly uncertain. The correction is about 3%, < 1%, 4% and < 2% for the 217 × 545, 353 × 545, 217 × 857 and 353 × 857 cross-power spectra, respectively. Although small, for consistency with the case  GHz and  GHz, we also apply the tSZ-related corrections to the measured CIB when  GHz and  GHz.

3.2. Dust model

Many studies, using mostly IRAS and COBE data, have revealed the strong correlation between the far-infrared dust emission and 21-cm integrated emission at high Galactic latitudes. In particular, Boulanger et al. (1996) studied this relation over the whole high Galactic latitude sky and reported a tight dust-Hi correlation for N_HI< 4.6 × 10²⁰ cm^-2. For higher column densities the dust emission systematically exceeds that expected by extrapolating the correlation. Examining specific high Galactic latitude regions, Arendt et al. (1998), Reach et al. (1998), Lagache et al. (1998 and 1999) found infrared excesses with respect to N_HI, with a threshold varying from 1.5 to 5.0 × 10²⁰ cm^-2. Planck Collaboration 2011XPlanck Collaboration XXIV (2011) presented the results from the comparison of Planck dust maps with GBT Hi observations in 14 fields covering more than 800 deg². They showed that the brighter fields in their sample, with an average Hi column density greater than about 2.5 × 10²⁰ cm^-2, show significant excess dust emission compared to the Hi column density. Regions of excess lie in organized structures that suggest the presence of hydrogen in molecular form. Because of this, we restrict our CIB analysis to the cleanest part of the sky, with mean N_HI< 2.5 × 10²⁰ cm^-2.

3.2.1. Constructing dust maps

As detailed in Planck Collaboration 2011XPlanck Collaboration XXIV (2011), Planck Collaboration 2011RPlanck Collaboration XVIII (2011) and in Appendix A, we constructed integrated Hi emission maps of the different Hi velocity components observed in each individual field: the local component, typical of high-latitude Hi emission, intermediate-velocity clouds (IVCs), and high-velocity clouds (HVCs), if present. To remove the cirrus contamination from HFI maps, we determined the far-IR to millimetre emission of the different Hi components. We assumed that the HFI maps, I_ν(x,y), at frequency ν can be represented by the following model ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}}{\mathit{\alpha }}_{\mathit{\nu }}^{\mathit{i}}{\mathit{N}}_{\mathrm{HI}}^{\mathit{i}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathit{\nu }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathit{,}$(13)where ${\mathit{N}}_{\mathrm{HI}}^{\mathit{i}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}$ is the column density of the ith Hi component, ${\mathit{\alpha }}_{\mathit{\nu }}^{\mathit{i}}$ is the far-IR to mm versus Hi correlation coefficient of component i at frequency ν and C_ν(x,y) is an offset. The correlation coefficients ${\mathit{\alpha }}_{\mathit{\nu }}^{\mathit{i}}$ (often called emissivities) were estimated using χ² minimization given the Hi, HFI and IRIS data, as well as the model (Eq. (13)). We removed from the HFI and IRIS maps the Hi velocity maps multiplied by the correlation coefficients. For EBHIS and GBT fields, we only considered one correlation coefficient per field and per frequency. The removal was done at the HFI and IRIS angular resolutions, even though the Hi map is of lower resolution (~10′). This is not a problem because cirrus, with a roughly k^-2.8 power-law power spectrum (Miville-Deschênes et al. 2007), has negligible power between the Hi (GBT and EBHIS) and HFI and IRIS angular resolutions, in comparison to the power in the CIB. The correlation of the dust emission with the different Hi velocity components and its variation from field to field is illustrated in Fig. 5 of Planck Collaboration 2011RPlanck Collaboration XVIII (2011).

For the GASS field, due to its large size the dust model needs to take into account variations of the dust emissivity across the field. We make use of an analysis of the dust-to-gas correlation over the southern Galactic cap (b< − 30°) using Hi data from the GASS southern sky survey (Kalberla et al. 2010). The Planck-HFI maps are linearly correlated with Hi column density over an area of 7500 deg² covering all of the southern sky (δ< 0°) at b< − 30° (17% of the sky). We use HFI maps corrected for the mean value of the CIB (Planck Collaboration 2013HPlanck Collaboration VIII 2014). The Planck maps and Hi emission at Galactic velocities are correlated over circular patches with 15° diameters, centred on a HEALPix grid with Nside = 32. The linear regression is iterated to identify and mask sky pixels that depart from the correlation. At microwave frequencies the correlation coefficients (α_ν) and offsets (C_ν) derived from this linear correlation analysis include a significant CMB contribution that comes from the chance correlation of the cosmic background with the Hi emission. This contribution is estimated by assuming that the SED of dust emission follows a modified blackbody spectrum for 100 ≤ ν ≤ 353 GHz. The fit is performed on the differences α_ν − α_{100 GHz} that are CMB-free for each sky area when expressed in units of K_CMB. This yields values of the correlations coefficients corrected for CMB, ${\mathit{\alpha }}_{\mathit{\nu }}^{\mathrm{c}}$. The detailed procedure is described in Planck Collaboration (2013). In this section we explain how the results of this study are used to build a model of the dust contribution to the sky emission.

To make the dust model in the GASS field, we start by building a map of the dust emission to Hi column density ratio, interpolating the values of ${\mathit{\alpha }}_{\mathit{\nu }}^{\mathrm{c}}$, corrected for the CMB, using HEALPix sky pixels with a Gaussian kernel. The 1σ width of this convolution kernel is equal to the pixel size 1.8° of the HEALPix grid for Nside = 32. To reduce the data noise at ν< 353 GHz we use the modified blackbody fits to α_ν − α_{100 GHz} and not the measured values of α_ν. This yields a set of six maps of the dust emission per unit Hi column density for all HFI frequencies from 100 to 857 GHz. We also build a set of Galactic offset maps from the offsets C_ν of the Planck-Hi correlation. These offsets comprise contributions from Galactic dust and the CMB. We subtract the CMB contribution assuming that the SED of the dust contribution to C_ν is the same as that of ${\mathit{\alpha }}_{\mathit{\nu }}^{\mathrm{c}}$ for each sky area. For each frequency, the dust model is the product of the dust emission per unit Hi column density times the Hi map, plus the Galactic offset map. The angular resolution of the model is that of the Hi map (16.2′).

As for the smaller EBHIS and GBT fields, the model in the GASS field only accounts for the emission of dust in Hi gas. Clouds with a significant fraction of molecular gas produce localized regions with positive residual emission. The histogram of residual emission at 857 GHz also shows a non-Gaussian extension towards negative values (as observed in the EBHIS field, see Appendix A.3). In the maps these pixels correspond to localized Hi clouds with no (or a weak) counterpart in the Planck map. In the GASS survey, these clouds are likely to be part of the Magellanic stream, with radial velocities within the range used to build the Galactic Hi map. For the analysis of the CIB we mask pixels with positive and negative residuals larger than 3σ. To be conservative this first mask is slightly enlarged, and apodized. It covers about 4400 deg² (Mask2 in Table 2). This mask still contains some regions with N_HI ≥ 2.5 × 10²⁰ cm^-2 and thus potential IR emission from molecular gas clouds. To measure the CIB power spectrum, all regions with N_HI ≥ 2.1 × 10²⁰ cm^-2 are further masked (Mask1). The final area is about 1900 deg².

3.2.2. Dust map uncertainties

The uncertainties estimated for the emissivities by the least-squares fit method are substantially underestimated, as they do not take into account systematic effects associated with the CIB and the CMB, nor the spatial variation of the emissivity inside a patch (or a field). We use the GASS field to estimate the error we have on the dust model. For this field, we have the large-scale variations (at ~15°) of the dust emissivities. We make the hypothesis that the measured variations extend to smaller angular scales, with a distribution on the sky given by a power spectrum with a slope − 2.8, similar to that of the dust emission. Within this assumption, we simulate multiple maps of the dust emissivity that all match the dispersion of the dust emissivity measured at 857  GHz. For each realization, we obtain a dust map at 857  GHz by multiplying the dust emissivity map by the GASS map of Galactic Hi. Dust maps at other frequencies are obtained by scaling the 857  GHz flux using the mean dust SED given. These simulations provide a good match to the Galactic residuals of the dust-Hi correlation characterized. We obtain simulated maps of the sky emission adding realizations of the HFI instrument noise, CIB and CMB to the dust maps. We perform on the simulated sky maps the same correlation analysis with Hi as that done on the Planck maps. For each simulation, we obtain values of the dust emissivity that we compare with the input emissivity map averaged over each sky patch. We find that there is no systematic difference between the values derived from the correlation analysis and the input values of the dust emissivities. The fractional error, i.e., the standard deviation of the difference between measured and the input values divided by the mean dust emissivity, is 13% of the mean dust emissivity at 857  GHz. This error increases slightly towards lower frequencies up to 16 and 21% for the 143 and 100  GHz channels.

We show in Fig. 3 the power spectrum of the recovered CIB, compared to the input CIB. We see that it is strongly biased by Galactic dust residuals in the first two multipole bins. These points have thus to be considered as upper limits. For the other points, we use the simulations to set the error bars linked to the Galactic dust removal. The observed linear correlation between the sigmas of dust residuals and N_HI is used to compute the power spectrum of Galactic dust residuals for each field, following ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{Field}}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{GASS}}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\left(\frac{\mathrm{⟨}{\mathit{N}}_{\mathrm{HI}}^{\mathrm{Field}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{N}}_{\mathrm{HI}}^{\mathrm{GASS}}\mathrm{⟩}}\right)}^{\mathrm{2}}\mathrm{·}$(14)

3.2.3. Galactic dust residuals

One of the main issues when using Hi column density as a dust tracer is the presence of emissivity variation on scales smaller that those probed by the correlation analysis, “dark gas” (see Planck Collaboration 2011SPlanck Collaboration XIX 2011; Planck Collaboration 2011XPlanck Collaboration XXIV 2011), and ionized gas (Lagache et al. 2000). The dust contribution of the dark gas becomes rapidly visible (not only at high frequencies), when N_HI ≳ 2.5 × 10²⁰ cm^-2. In addition to the simulations discussed in the previous section, we investigated the contribution of dust residuals by computing the CIB power spectra on the GASS field, using Mask1 and Mask2. Mask1 is very conservative and is our nominal mask for CIB analysis. Mask2 contains higher column density regions, so it is not suitable for CIB-only analysis, but is useful, for example, when cross-correlating the CIB with other large-scale structure tracers (when dust contamination is less of a problem, e.g., Planck Collaboration 2013RPlanck Collaboration XVIII 2014) and for the bispectrum measurement at low frequencies. As expected, the power spectra of residual maps computed using Mask2 have an excess compared to the CIB. Whatever the frequency, this excess represents about 5–10% of the cirrus power spectrum at low ℓ. This leads to a CIB from Mask2 that is overestimated by a factor of 1.5 at 857  GHz and 1.2 at 217  GHz, for ℓ ≤ 200, compared to the CIB from Mask1. Taking such a residual at the level of 5–10% of the dust power spectrum (as observed between Mask1 and Mask2) would strongly affect only the first two bins, at ℓ = 53 and 114. This is consistent with the simulation analysis (see Fig. 3). We will thus not use those bins when searching for the best CIB model, and we consider them only as upper limits.

3.3. CMB and Galactic dust cleaned maps

We show in Fig. 4, one example of the cleaning process, from the frequency to the CIB maps, for the SPC5 field. The bottom row shows the residual CIB maps, smoothed to 10′. In Fig. 5 are shown the CIB maps in a part of the GASS field, at 16.2′ resolution. We see from both figures that common structures, corresponding to CIB anisotropies, are clearly visible. As previously noticed (Planck Collaboration 2011RPlanck Collaboration XVIII 2011), the three intermediate frequencies (353, 545 and 857  GHz) show highly correlated structures. On the contrary, the 3000  GHz data on one hand, and the 217  GHz on the other, reveal a decoherence, which can be attributed to the redshift distribution of the CIB anisotropies. We will come back to this decoherence by measuring the correlation coefficients in Sect. 6.2.

There are two frequencies where extracting the CIB from the frequency maps is particularly challenging.

• At 3000 GHz, the cirrus contamination is the highest and we observe more spatial variations of the dust–Hi emissivity which are difficult to manage. This is due to the contamination at 3000 GHz data by a hotter dust component, which may be linked to the so-called “very small grains” in some interstellar filaments. Consequently, we have not tried to extract the CIB in the GASS field at 3000 GHz, and the Bootes field has been discarded for the CIB analysis due to interstellar dust residuals (due to one IVC component, which has both a high emissivity and a hot spectrum; P. Martin, priv. comm.). Moreover, the EBHIS field is right in the middle of the missing IRAS observation. It is thus also not used for the CIB analysis. In the end, the total area used to compute the CIB power spectrum at 3000 GHz is 183 deg².

• At 143 GHz, we expect the correlated CIB anisotropies in brightness to be about 1–2% of the CMB anisotropies for multipoles ℓ ~ 100–1000 (while it is about 3–15% at 217GHz). The removal of the CMB has thus to be extremely accurate. Moreover, the expected CIB is lower than the instrument noise. Constraints can only be obtained on the large GASS field that is more immune to noise.We see from Fig. 5 that the 143GHz CIB map shows some structures that are correlated with the higher CIB frequency maps. We discuss in Sect. 6.1 our attempt to obtain some constraints at this frequency.

We used two different approaches to measure the CIB power spectra according to the size of the fields: (i) for the EBHIS and GBT flat-sky fields, we used an updated version of POKER (Ponthieu et al. 2011), and we computed the error bars using Monte Carlo simulations (Sect. 4.1); (ii) for the GASS field, we used Xspect (Tristram et al. 2005), a method that was first developed for the Archeops experiment to obtain estimates of the angular power spectrum of the CMB temperature anisotropies, including analytical error bars (Sect. 4.2). Having two completely different pipelines and a many fields with various dust and CMB contaminations (and noise contributions), is extremely valuable, as it allows us to test the robustness of our approach.

4. Angular power spectrum and bispectrum

4.1. Cross-correlation pipeline for the EBHIS and GBT fields

To determine the cross-correlation between the CIB observed at two frequencies, ν₁ and ν₂, we used a modified version of POKER (Ponthieu et al. 2011). POKER is an algorithm that determines the power spectrum of a map, corrects for mask aliasing and computes the covariance between each power spectrum bin via a Monte Carlo approach. It has been used to measure the power spectrum of the CIB as observed by Planck-HFI (Planck Collaboration 2011RPlanck Collaboration XVIII 2011) and is described in detail therein. In the following, we will call auto-power spectrum ${\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}}$ the usual angular power spectrum and cross-power spectrum ${\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}$ its generalization to two different frequencies defined as $\frac{\mathrm{1}}{\mathrm{2}}\mathrm{⟨}{\mathit{a}}_{\mathit{\ell m}}^{{\mathit{\nu }}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{{\mathit{\nu }}_{\mathrm{2}}\mathrm{\ast }}\mathrm{+}{\mathit{a}}_{\mathit{\ell m}}^{{\mathit{\nu }}_{\mathrm{1}}\mathrm{\ast }}{\mathit{a}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{{\mathit{\nu }}_{\mathrm{2}}}\mathrm{⟩}\mathrm{=}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{,}$(15)where ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{\nu }}$ are the Fourier coefficients of the CIB anisotropy at the observation frequency ν. We take a common mask for both frequencies and it is straightforward to generalize the algebra presented in Ponthieu et al. (2011) and Planck Collaboration 2011RPlanck Collaboration XVIII (2011) to obtain an unbiased estimate of the cross-power spectrum: ${\mathit{C}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}\mathrm{\simeq }\sum _{{\mathit{b}}^{\mathrm{\prime }}}{\mathit{M}}_{\mathit{b}{\mathit{b}}^{\mathrm{\prime }}}^{-1}\mathit{P̂}\begin{array}{c}{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\\ {\mathit{b}}^{\mathrm{\prime }}\end{array}\mathrm{-}{\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{instr}}\mathrm{-}{\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{res}}\mathit{.}$(16)Here $\mathit{P̂}\begin{array}{c}{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\\ \mathit{b}\end{array}$ denotes the binned pseudo-cross-power spectrum of the maps at frequencies ν₁ and ν₂, and M_bb′ is the so-called mode-mixing matrix that is described in detail in Eq. (24) of Planck Collaboration 2011RPlanck Collaboration XVIII (2011) and its appendix. We recall that it includes the mode-coupling effects induced by the mask, the smoothing by the instrumental beam, and the map pixelization (see below). The two noise terms, ${\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{instr}}$ and ${\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{res}}$, refer, respectively, to the instrument noise and to the contribution of what can be considered as random components at the map level, such as the CMB residuals of our component separation and the noise of the 100 GHz map that propagates to our maps (see Sect. 3.1). Of course, Eq. (16) applied to ν₁ = ν₂ gives the auto-power spectrum of the anisotropy, and we computed the two auto-spectra at ν₁ and ν₂ at the same time as the cross-power spectrum.

We built very similar simulation and analysis pipelines to those of Planck Collaboration 2011RPlanck Collaboration XVIII (2011) to obtain our measures and their associated error bars. We here briefly recall the main steps of these pipelines and highlight the modifications introduced by the generalization from auto-power spectra estimation to that of cross-power spectra.

Analysis pipeline.

• 1.

  In order to combine our measurements from different fields, wedefine a common multipole binning. Above ℓ ~ 200, we choose a logarithmic binning, Δℓ/ℓ = 0.3, while below ℓ ~ 200, we respect the generic criterion that bins of multipoles should be larger than twice the multipole corresponding to the largest angular scale contained in the field.

• We define a common weight map W^ν₁ν₂ that masks out bright point sources found at both frequencies. Together with the mask, three different transfer functions must be accounted for in the computation of M_bb′ (Eq. (24) of Planck Collaboration 2011RPlanck Collaboration XVIII 2011): (i) the instrument beam transfer function that depends on the exact beam shape and the scan pattern on the observed field, although since the variation of beam transfer function is less than 1% between our fields, we take the same one for all of them (see Sect. 2.1); (ii) the HEALPix pixel window function; and (iii) the reprojection from HEALPix maps to the flat-sky maps used by POKER. The first transfer function comes from a dedicated analysis (Planck Collaboration 2013GPlanck Collaboration VII 2014). The second one is provided by the HEALPix library. The third one is computed via Monte Carlo simulations: we simulate full sky HEALPix maps of diffuse emission with a typical CIB power spectrum, reproject them on our observed patches and compute the angular power spectra; the ratio between the measured and the input power spectrum gives the transfer function. This ratio is the same for all fields (all in SFL projection with 3.5′ pixels) within statistical error bars, so we compute the average and apply it to all fields.

  Table 4

  Example of a power spectrum averaged for the flat fields.

• 3.

  An estimate of the noise auto-power spectrum at ν₁ and ν₂ is obtained from jack-knife maps. Indeed, at each frequency, two maps can be built using only the first (respectively, the second) half of each Planck observation ring. The difference between these maps is dominated by instrumental noise. Applying POKER to these difference maps gives an estimate of the instrument noise auto-power spectrum at each frequency ${\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathit{i}}\mathit{,}\hspace{0.17em}\mathrm{instr}}$. The noise contribution to the auto-power spectrum are negligible at high frequencies.

• 4.

  ${\mathit{N}}_{\mathit{b}}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,}\hspace{0.17em}\mathrm{res}}$ is the contribution of CMB residuals both from the component separation and from the propagation of noise in the 100 GHz CMB map. We have estimates of each component of this residual power spectrum in ℓ space on the sphere (see. Sect. 3.1) and combine them into our measurement bins.

• 5.

  We now apply Eq. (16) to our data and determine their auto- and cross-power spectra.

Simulation pipeline.

The simulation pipeline is essential to provide the error bars on our estimates coming from the analysis pipeline. We created 100 simulations of our data maps and computed their auto- and cross-power spectra. The dispersion of these spectra gives the complete covariance matrices of our binned auto- and cross-power spectra and their associated error bars. For each realization we follow these steps.

• 1.

  The measured auto- and cross-power spectra are used as inputsto simulate CIB anisotropy maps at each frequency, with theappropriate correlated component. To do so, we generate randomGaussian amplitudes x_ℓ and y_ℓ in Fourier space, such that: $\begin{array}{ccc}{\mathit{a}}_{{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}}& \mathrm{=}& {\mathit{x}}_{{\ell }}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{;}\\ {\mathit{a}}_{{\ell }}^{{\mathit{\nu }}_{\mathrm{2}}}& \mathrm{=}& {\mathit{x}}_{{\ell }}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}\mathit{/}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{+}{\mathit{y}}_{{\ell }}{\left\{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{2}}{\mathit{\nu }}_{\mathrm{2}}}\mathrm{-}{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{/}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}_{\mathrm{1}}}\right\}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$

• 2.

  For each frequency ν_i, we simulate a noise map with the appropriate power spectrum $\mathit{N̂}\begin{array}{c}{\mathit{\nu }}_{\mathit{i}}\mathit{,}\hspace{0.17em}\mathrm{instr}\\ \mathit{\ell }\end{array}$ and add it to the simulated signal map.

• 3.

  We add the simulations of CMB and 100  GHz noise residuals (Eqs. (1) and (2)), and of CMB residuals induced by relative calibration errors, so that we have a final pair of maps at ν₁ and ν₂ that are a faithful representation of our data.

Error estimation.

Statistical uncertainties due to instrument noise and component separation residuals are derived using the simulation pipeline. Systematic uncertainties include mask aliasing effects, imperfect subtraction of foreground templates, and beam and projection transfer function errors. We studied some fields in common with Planck Collaboration 2011RPlanck Collaboration XVIII (2011), with very similar masks, and showed that the correction of mask effects by POKER leads to no more than 2% uncertainty on the result.

The transfer function due to the projection of spherical HEALPix maps to our square patch is determined via Monte Carlo (see step 2 of the description of the analysis pipeline). This process provides estimates F_b of the transfer function in our measurement bins. However, we need an estimate of the transfer function F_ℓ for all ℓ modes, to include it in the derivation of M_bb′. We therefore construct a smooth interpolation of our measures. The projection transfer function plays the same role as an extra instrumental beam and enters the derivation of M_bb′. Roughly speaking, it damps the signal at high ℓ and correcting for this effect corresponds to dividing the measured power spectrum by the damping function. The statistical uncertainty on the determination of F_ℓ is smaller than 1%, and hence leads to the same uncertainty on the C_ℓ.

Between the first analysis of CIB anisotropies with Planck-HFI and this work, much progress has been achieved on the beam measurements (Planck Collaboration 2013GPlanck Collaboration VII 2014). Based on this progress, we have better estimates of the beam shapes and better assessments of the uncertainties on these beam transfer functions, both for auto-spectra and for cross-spectra between different frequency bands. We used the eigenmodes to compute the beam transfer function uncertainties, which are all smaller than 0.55% in our angular range. On the contrary, auto- and cross-spectra involving IRIS suffer from a much larger beam uncertainty (0.5′ for a FWHM of 4.3′).

The systematic uncertainties on the contribution of Galactic dust residuals (see Sect. 3.2.2) are added linearly to the statistical error bars. An example of the error budget is given in Table 4.

4.2. Cross-correlation pipeline for the GASS field

For the large GASS field, we used another strategy:

• Due to the size of the field, which violates the flat-skyapproximation used for the EBHIS and GBT fields, we computethe angular auto- and cross-power spectra on theHEALPix maps, using the Xspect algorithm(Tristram et al. 2005). We apply theGalactic dust mask discussed in Sect. 3.2.1. Thepower spectra are computed using maps at the HFI angularresolution (N_side = 2048).

• As the angular resolution of the Hi map is 16′, we use a hybrid method to remove the Galactic dust. On large angular scales (ℓ ≤ 590), we remove the cirrus from the maps using Hi. On small angular scales (ℓ > 590), we remove an estimate of the dust power spectrum. This estimate comes from the dust model that is fit on large angular scales (120 ≤ ℓ ≤ 590) by a power law, following Miville-Deschênes et al. (2007), and then extrapolated to small angular scales. The power-law fit has been shown to be valid for the whole range of angular scales covered by our measurements (e.g., Miville-Deschênes et al. 2010).

• As for the flat-sky fields, the CMB is removed as described in Sect. 3.1.

• Statistical error bars are not computed using simulations but using analytical formulations (see below).

Auto- and cross-power spectra.

We used the maps built from the first and second halves of each pointing period. As described in Tristram et al. (2005), the spherical harmonic coefficients from the cut-sky maps are corrected from the mode-coupling introduced by the mask, as well as the beam smoothing effect in the harmonic domain. The cross-power spectra are unbiased estimates of the angular power spectrum, avoiding any correction for the instrument noise, contrary to our measurements done in Sect. 4.1. We computed the spectra using the same multipole binning as that of the flat fields. We checked that the bin to bin correlation was smaller than 1%.

We show in Fig. 6 an example of a cross-power spectrum (545 × 353) computed on maps for which the CMB has been removed. We also show the power spectrum of the dust model, fitted by a broken power law. The slope of the dust power spectrum appears to flatten at ℓ ~ 110, with a slope in the range − 2.1 to − 1.8 for ℓ < 110, and − 2.8 to − 2.7 for ℓ > 110, depending on the frequency. We compare the CIB power spectrum obtained by removing the dust at the map level to that obtained by subtracting the fit of the dust model power spectrum. For ℓ ≥ 200 we have excellent agreement. We also see that for ℓ ≥ 700 the dust removal has a negligible impact. At low ℓ, where the dust correction is important, we chose to use the dust removal on the map, since it leads to a lower variance on the residual power spectrum, as shown in Pénin et al. (2012b). Indeed, the spatial subtraction removes each moment of the statistics, whereas the subtraction of the power spectrum only removes the first two moments. At low ℓ, errors on the fit are also quite large. We arbitrary decide to take the transition between the dust removal on the map and on the power spectrum at ℓ = 510. The exact choice of the multipole for the transition has no consequence on the resulting CIB power spectrum.

Error bars.

We quadratically combine the following error terms to obtain the final uncertainty on the GASS CIB power spectra:

• The error bars on the spectrum computed analytically asdescribed in Tristram et al. (2005),from the auto-power and cross-power spectra of the two maps.They include both the sampling variance (which dominates atlarge scales) and the instrumental noise (which dominates atsmall scales).

• The errors linked to the CMB removal. Errors on the extra instrumental noise that have been introduced by the CMB removal (see Eq. (1)) are computed using the errors on the noise measurements at 100 GHz. For the error on the CMB contribution that is left close to the angular resolution of the 143 and, 217 and 353 GHz channels (see Eq. (2)), we use the theoretical cosmic variance estimate on the determination of the CMB power spectrum. Finally, we add the errors that come from the relative photometric calibration in the CMB removal.

As detailed in Sect. 3.2.2, the uncertainty on the dust model is derived from simulations. It is added linearly to the statistical error bars. For ℓ ≥ 510, we also add the error on the dust model fit. Beam errors are computed using the eigenmodes. They are the same as those detailed in Sects. 2.1 and 4.1.

4.3. Bispectrum pipeline for the GASS field

The bispectrum b_ℓ1ℓ2ℓ3 is the 3-point correlation function in harmonic space: $\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{=}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{G}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}\mathit{,}$(17)with ${\mathit{G}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}$ the Gaunt integral (Spergel & Goldberg 1999). It is a lowest-order indicator of the non-Gaussianity of the field.

The maps used are the same as for the power spectrum analysis, but degraded to N_side = 512. As the S/N for bispectra is quite low compared to that of power spectra, bispectra have to be measured on the largest possible clean area of the sky. Here, we measure the bispectrum on GASS using Mask2. We apply the binned bispectrum estimator described in Lacasa et al. (2012) and used for the Planck tSZ map analysis (Planck Collaboration 2013UPlanck Collaboration XXI 2014) and non-Gaussianity constraints (Planck Collaboration 2013XPlanck Collaboration XXIV 2014). A large bin size Δℓ = 128 has been adopted to minimize multipole correlations due to the mask. We used a multipole range ℓ_min = 129 to ℓ_max = 896, leaving six multipole bins and 43 bispectrum configurations (ℓ₁,ℓ₂,ℓ₃) (accounting for permutation invariance and the triangular condition). We only considered auto-bispectra for simplicity, i.e., at a single frequency (see Table D.3).

For each frequency, we computed the auto-bispectra of the two maps built using the two half-pointing period rings, and average these bispectra for a raw estimate. The raw estimate has been then debiased from mask and beam effects using simulations. Specifically we generated simulations with a high level of non-Gaussianity and a bispectrum corresponding to the CIB prescription from Lacasa et al. (2012), computed the ratio of the bispectrum of the masked map to the full-sky bispectrum, and finally averaged this ratio over simulations, finding a quick convergence especially at high multipoles. We found this ratio to be very close to f_skyb_ℓ1(ν) b_ℓ2(ν) b_ℓ3(ν) (with b_ℓ(ν) the beam window function), showing that the multipole correlation is indeed negligible for this bin size and bispectrum. We checked on the two half maps that the Planck noise is close to Gaussian, hence it does not bias our bispectrum estimates; it however increases their variance.

The error estimates are the sum in quadrature of:

• cosmic variance computed with analytical formulae andincluding the noise;

• dust residuals from (the absolute value of) the bispectrum of the dust model, scaled to the residual dust amplitude found in Sect. 3.2.2.

Note that beam errors are completely negligible in the range of ℓ considered. The full-sky cosmic variance of the bispectrum is composed of four terms: $\mathrm{Cov}\mathrm{(}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\hspace{0.17em}\mathit{,}\hspace{0.17em}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}\mathrm{)}\mathrm{=}{\mathrm{C}}_{\mathrm{2}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{2}}\mathrm{+}{\mathrm{C}}_{\mathrm{3}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{3}}\mathrm{+}{\mathrm{C}}_{\mathrm{2}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{4}}\mathrm{+}{\mathrm{C}}_{\mathrm{6}}\mathit{,}$(18)with (see Lacasa et al. 2012) ${\mathrm{C}}_{\mathrm{2}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{2}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{2}}\mathrm{=}\frac{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}\hspace{0.17em}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}\hspace{0.17em}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}}{\frac{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}{\mathrm{)}}_{\mathrm{123}}}{\mathrm{4}\mathit{\pi }}{\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& {\mathit{\ell }}_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)}^{\mathrm{2}}}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{\Delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathit{.}$(19)Here ${\mathrm{\Delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{6}& \mathrm{equilateral}\mathrm{triangles}\\ \mathrm{2}& \mathrm{isosceles}\mathrm{triangles}\\ \mathrm{1}& \mathrm{general}\mathrm{triangles}\end{array}$(20)and when the bispectrum is factorizable $\begin{array}{ccc}{\mathrm{C}}_{\mathrm{3}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{3}}\mathrm{=}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\hspace{0.17em}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}& \mathrm{\times }& (\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}}\\ & & \mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}}\\ & & \mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}}\mathrm{+}\frac{{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}}}{\mathrm{2}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}})\end{array}$(21)C_{2 × 4} and C₆ involve, respectively, the map trispectrum and 6-point function. C_{2 × 2 × 2} is the cosmic variance in the weak non-Gaussianity limit, as considered, e.g., by Crawford et al. (2014) in their tSZ and CIB non-Gaussianity study; it produces a diagonal covariance matrix, since it does not couple the different configurations. C_{3 × 3} is the bispectrum correction to the covariance matrix and couples bispectrum estimates in different configurations. In the cosmic variance limited case, we found that including only C_{2 × 2 × 2} would noticeably overestimate the detection significance. However when considering all error sources, C_{3 × 3} has a small impact on the S/N; hence we include it in the following analysis but neglect higher order corrections C_{2 × 4} and C₆.

We use the f_sky approximation, that is Cov = Cov_{full − sky}/f_sky, as the first multipoles were discarded and we saw no mode coupling with our large bin size. The power spectrum used for C_{2 × 2 × 2} is the map auto-power spectrum, including the noise as necessary, debiased from the mask and beam effects.

As described in Sect. 3.1.2, the residual maps are contaminated by tSZ contribution at 100  GHz leaking through the CMB cleaning process. While this contamination is negligible compared to the CIB signal at 353  GHz and above, it is important at 217  GHz, where the contamination can be as large as 40%, depending on the configuration and multipole. We derived the SZ bispectrum by measuring the bispectrum of the map of the Planck SZ cluster catalogue (both clusters and candidate clusters from Planck Collaboration 2013ZDPlanck Collaboration XXX (2014) are considered). An 80% error on the amplitude of the tSZ correction is added to the covariance measurement described previously. This error is based on the relative difference between the bispectrum of the Planck SZ catalogue of confirmed clusters (Planck Collaboration 2013ZDPlanck Collaboration XXX 2014), the bispectrum of the Planck estimated tSZ map, and the bispectrum of the Planck FFP6 SZ simulations (see the bispectrum analysis in Planck Collaboration 2013UPlanck Collaboration XXI 2014). This error estimate is conservative, since it takes the maximum observed difference for all multipoles and configurations.

4.4. CIB power spectrum

Figure 7 presents a summary of all the measured auto- and cross-power spectra on residual maps. Following the same covariance studies as in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), we combine our cross-power spectrum estimates on individual fields f for each bin b into an average cross-power spectrum, using inverse variance weights, ${\mathit{C}}_{\mathit{b}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}\frac{{\sum }_{\mathit{f}}{\mathit{W}}_{\mathit{b}}^{\mathit{f}}{\mathit{P}}_{\mathit{b}}^{\mathit{f,\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}}{{\sum }_{\mathit{f}}{\mathit{W}}_{\mathit{b}}^{\mathit{f}}}\mathit{,}$(22)with ${\mathit{W}}_{\mathit{b}}^{\mathit{f}}\mathrm{=}\mathrm{1}\mathit{/}{\mathit{\sigma }}^{\mathrm{2}}\left({\mathit{P}}_{\mathit{b}}^{\mathit{f,\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\right)$. These weights are estimated in the following way:

• Step 1: only statistical errors are used to compute a weightedaverage of the power spectra and their associated error bars onsmall fields.

• Step 2: the projection error is added linearly to the error bar of step 1.

• Step 3: for frequencies where observations on GASS are available, we average the GASS power spectrum and the small fields power spectrum. We here use inverse statistical variance weights for GASS, and the inverse variance derived from step 2 for the small fields. This gives the final power spectrum estimate and a pseudo-statistical error bar.

• Step 4: we compute the error due to beam uncertainties using the average power spectrum. We add this error linearly to the error derived on Step 3.

• Step 5: the bias induced by Galactic dust residuals is linearly added to the error derived on Step 4 to obtain the final total error bar.

The resulting power spectra and their errors are given in Table D.1. For the 217 × 217 auto-power spectrum, only the measurement on the GASS field is considered, as the measurement in the flat-sky fields is noise dominated. Note also that the last bin for all measurements involving the subtraction of the 100  GHz CMB template is ℓ = 1956, due to the angular resolution of Planck-HFI at 100  GHz. To obtain the CIB power spectra, the estimates obtained from the CMB- and dust-free maps (Table D.1) have to be corrected for SZ contamination (Eqs. (5) and (9)), and for the spurious CIB contamination (Eq. (4)) induced by our CMB template. This last contaminant is computed using our best-fit model described in Sect. 5.5. Moreover, following Sects. 3.2.2 and 3.2.3, the first two bins at multipoles ℓ = 53 and 114 have to be considered as upper limits. CIB power spectrum values are given in Table D.2. The errors contain all the terms: statistical uncertainty; beam and projection uncertainty; cirrus bias; and errors from the SZ correction. CIB power spectra, comparing the measurements with the model are presented later (Fig. 12). Comparison with previous recent measurements are shown and discussed in Sect. 6.3.

The power spectra from the flat-sky and GASS fields have been obtained using two independent pipelines. The fields have different Galactic dust and point sources contamination, as well as instrument noise levels, and between GASS and the flat-sky fields, different pixelization and projection. Comparing combined power spectra obtained on all flat-sky fields to that obtained on GASS is thus a powerful consistency check on our determination of power spectra and error bars. We show in Fig. 8 this comparison for an arbitrary set of frequencies. The power spectra are always compatible within 1σ. Of course, due to the much larger area of the GASS field, the CIB measured in GASS has much smaller cosmic variance errors.

4.5. CIB bispectrum

We measure the bispectrum only at 217, 353, and 545  GHz. The bispectrum at 143  GHz is noise-dominated and is moreover highly contaminated by tSZ and extragalactic radio point sources. At 857  GHz, as we are using Mask2, Galactic dust residuals contaminate the bispectrum in most configurations. In particular, the residuals produce a rising bispectrum at high multipoles. In the following, we are thus not considering the 857  GHz frequency for the analysis. At 217, 353, and 545  GHz, the bispectrum is measured in 38, 40 and 36 configurations, respectively. We show in Fig. 9 the measured CIB bispectrum at 353  GHz for some particular configurations, namely equilateral (ℓ,ℓ,ℓ), orthogonal isosceles $\mathrm{\left(}\mathit{\ell ,\ell ,}\sqrt{\mathrm{2}}\mathit{\ell }\mathrm{\right)}$, flat isosceles (ℓ,ℓ,2ℓ), and squeezed (ℓ_min,ℓ,ℓ). The bispectrum decreases with scale and exhibits a peak in the squeezed configurations, as predicted by Lacasa et al. (2012).

In Table 5, we give the significance of the detection for the three frequencies used, with either the average significance per configuration or the significance of the total bispectrum, when accounting for the whole covariance matrix. The bispectrum is significantly detected at each frequency individually. Moreover, these measurements represent the first detection of the CIB bispectrum per configuration, permitting us to probe the scale and configuration dependence of the bispectrum, as well as its frequency behaviour. The bispectrum values are given in Table D.3.

5. Interpreting CIB power spectrum measurements

Once the CMB and Galactic dust have been removed, there are three astrophysical contributors to the power spectrum at the HFI frequencies: two from dusty star-forming galaxies (with both shot noise, ${\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$, and clustering, ${\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{clust}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$, components); and one from radio galaxies (with only a shot-noise component, ${\mathit{C}}_{\mathrm{r}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$, the clustering of radio sources being negligible, e.g., Hall et al. 2010). The measured CIB power spectrum ${\mathit{C}}_{\mathrm{measured}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ is thus ${\mathit{C}}_{\mathrm{measured}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{clust}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{+}{\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{+}{\mathit{C}}_{\mathrm{r}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathit{.}$(23)In this section, we first discuss the shot-noise contributions. We then describe how we can model ${\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{clust}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$. Two approaches are considered. The first one is the simplest, and uses only the large-scale CIB measurements to fit for a linear model. The second one is based on the halo-model formalism. Our main goal is to use CIB anisotropies to measure the SFRD and effective redshift evolution of the bias.

5.1. Shot noise from dusty star forming and radio galaxies

The shot noise arises from sampling of a background composed of a finite number of sources, and as such is decoupled from the correlated term. The angular resolution of the HFI instrument is not high enough to measure the shot-noise levels. As demonstrated in Planck Collaboration 2011RPlanck Collaboration XVIII (2011), the non-linear contribution to the power spectrum is degenerate with the shot-noise level (on the scales probed by Planck). Since our data by themselves are not sufficient to explore this degeneracy, we need to rely on a model to compute the shot noise.

5.1.1. Auto-power spectrum

The shot-noise level at frequency ν can be easily computed using monochromatic galaxy number counts. Let us consider a flux interval [S_k,S_k + ΔS_k]. The number of sources per unit solid angle, n_k in this flux interval is ${\mathit{n}}_{\mathit{k}}\mathrm{=}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{\Delta }\mathit{S}k$(24)and the variance is, ${\mathit{\sigma }}_{{\mathit{B}}_{\mathit{k}}}^{\mathrm{2}}\mathrm{=}{\mathit{n}}_{\mathit{k}}{\mathit{S}}_{\mathit{k}}^{\mathrm{2}}\mathit{.}$(25)Summing all flux intervals gives the variance on the total contribution to the CIB: ${\mathit{\sigma }}_{\mathit{B}}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{k}}{\mathit{n}}_{\mathit{k}}{\mathit{S}}_{\mathit{k}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}{\mathit{S}}_{\mathit{k}}^{\mathrm{2}}\mathrm{\Delta }{\mathit{S}}_{\mathit{k}}\mathit{.}$(26)When we take the limit of ΔS_k tending to zero, this sum becomes the integral ${\mathit{\sigma }}_{\mathit{B}}^{\mathrm{2}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{S}}^{\mathrm{c}}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}{\mathit{S}}^{\mathrm{2}}\mathrm{d}\mathit{S}\mathit{.}$(27)Here S^c is the flux cut above which bright sources are detected and can be removed. This cut is mandatory, since the integral does not converge in the Euclidian regime, $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{\propto }{\mathit{S}}^{-2.5}$, which is the case at bright fluxes for star-forming galaxies.

5.1.2. Cross-power spectrum

We now consider two frequencies ν and ν′. The number of sources n_kl in the flux density and redshift intervals [S_k,S_k + ΔS_k] and [z_l,z_l + Δz_l], is ${\mathit{n}}_{\mathit{kl}}\mathrm{=}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}\mathrm{d}\mathit{z}}\mathrm{\Delta }{\mathit{S}}_{\mathit{k}}\mathrm{\Delta }{\mathit{z}}_{\mathit{l}}\mathit{.}$(28)Considering a small redshift interval, the covariance between the two frequencies can be approximated as: ${\mathit{\sigma }}_{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}\mathit{,kl}}\mathrm{=}{\mathit{n}}_{\mathit{kl}}{\mathit{S}}_{\mathit{\nu ,kl}}{\mathit{S}}_{{\mathit{\nu }}^{\mathrm{\prime }}\mathit{,kl}}\mathrm{=}{\mathit{n}}_{\mathit{kl}}{\mathit{S}}_{\mathit{\nu ,kl}}^{\mathrm{2}}{\mathit{R}}_{\mathit{\nu }{\mathit{\nu }}^{\mathrm{\prime }}\mathit{,kl}}\mathit{,}$(29)where R_νν′,kl is the mean colour for the considered galaxy population in the considered flux density and redshift interval. Using a mean colour per flux density and redshift interval is not a strong assumption as long as ΔS_k and Δz_l are small. Summing over flux densities, redshifts and the galaxy population gives ${\mathit{\sigma }}_{\mathit{\nu }{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}\sum _{\mathrm{pop}}\sum _{\mathit{k}}\sum _{\mathit{l}}{\mathit{n}}_{\mathrm{pop}\mathit{,kl}}{\mathit{S}}_{\mathit{\nu ,kl}}^{\mathrm{2}}{\mathit{R}}_{\mathit{\nu }{\mathit{\nu }}^{\mathrm{\prime }}\mathit{,}\mathrm{pop}\mathit{,kl}}\mathit{,}$(30)with the integral limit $\begin{array}{ccc}{\mathit{\sigma }}_{\mathit{\nu }{\mathit{\nu }}^{\mathrm{\prime }}}& \mathrm{=}& \sum _{\mathrm{pop}}{\mathrm{\int }}_{{\mathit{S}}_{\mathit{\nu }}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathrm{\int }}_{\mathit{z}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\mathit{H}\mathrm{\left(}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathit{<}\hspace{0.17em}{\mathit{S}}_{\mathit{\nu }}^{\mathrm{c}}\mathit{,}{\mathit{R}}_{{\mathit{\nu }}^{\mathrm{\prime }}\mathit{\nu ,}\mathrm{pop}}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathit{<}\hspace{0.17em}{\mathit{S}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{c}}\mathrm{\right)}\\ & & \end{array}$(31)Here H(P₁,P₂) is equal to 1 when P₁ and P₂ are both true, and 0 otherwise, and ${\mathit{S}}_{\mathit{\nu }}^{\mathrm{c}}$ and ${\mathit{S}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{c}}$ are the flux cuts in the frequency bands ν and ν′. They are given in Table 1.

We use this formalism to compute the radio galaxy shot noise. For the star-forming dusty galaxy shot noise, we rely on the formalism detailed in Béthermin et al. 2013 (their Appendix B).

5.1.3. Shot-noise values

We use the Béthermin et al. (2012a) model to compute the star-forming dusty galaxy shot noise, ${\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ (Eq. (23)). The model is in rather good agreement with the number counts measured by Spitzer and Herschel (e.g., Glenn et al. 2010). It also gives a reasonable CIB redshift-distribution, which is important for the cross-spectra. Since this model is based on observations that have typical calibration uncertainties of <8%, the estimations of the shot-noise levels (being proportional to the square of calibration factor) cannot be accurate to more than about 16%. Uncertainties in the flux cuts induce small variations in the shot noise (less than 3% at all frequencies). We take 20% as the shot-noise 1σ uncertainty.

For extragalactic radio sources, we use the Tucci et al. (2011) model (more specifically, the one referred as “C2Ex” in the paper) to compute ${\mathit{C}}_{\mathrm{r}\mathit{,}\mathrm{shot}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ (Eq. (23)). The predictions for high-frequency number counts are based on a statistical extrapolation of flux densities of radio sources from low-frequency data (1–5 GHz). In particular, this model considers physically based recipes to describe the complex spectral behaviour of blazars, which dominate the mm-wave counts at bright flux densities. It is able to give a good fit to all bright extragalactic radio source data available so far: number counts up to 600 GHz; and spectral index distributions up to at least 200–300  GHz (see Tucci et al. 2011; Planck Collaboration Int GPlanck Collaboration Int. VII 2013; López-Caniego et al. 2013). As for the dusty galaxies, we consider an error of 20% on the shot-noise computation from the model. But unlike for dusty galaxies, the shot noise for the radio population depends strongly on the flux cut (Planck Collaboration 2011RPlanck Collaboration XVIII 2011). Accordingly, we add to the 20% mentioned above, a shot-noise error deriving from the shot-noise variations as we change the flux cut, considering the flux cut errors given in Table 1.

The shot-noise levels for the HFI flux cuts given in Sect. 2.4, are listed in Tables 6 and 7.

5.2. Basics of CIB correlated anisotropy modelling

The angular power spectrum of CIB correlated anisotropies is defined as: $\begin{array}{ccc}\mathrm{⟨}\mathit{\delta }{\mathit{I}}_{\mathit{\ell m}}^{\mathit{\nu }}\mathit{\delta }{{\mathit{I}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{{\mathit{\nu }}^{\mathrm{\prime }}}}^{\mathrm{⟩}}\mathrm{=}{\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{clust}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{,}& & \end{array}$(32)where ν and ν′ denote the observing frequencies and I^ν,ν′ the measured intensity at those frequencies. In a flat universe, the intensity is related to the comoving emissivity j via $\begin{array}{ccc}{\mathit{I}}^{\mathit{\nu }}& \mathrm{=}& \int \mathrm{d}\mathit{z}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\mathit{aj}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{\int }\mathrm{d}\mathit{z}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\mathit{a}\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\left(\mathrm{1}\mathrm{+}\frac{\mathit{\delta j}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}}{\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}}\right)\mathit{,}\end{array}$(33)where χ(z) is the comoving distance to redshift z, and a = 1 / (1 + z) is the scale factor. Combining Eqs. (32) and (33) and using the Limber approximation, we obtain $\begin{array}{ccc}{\mathit{C}}_{\mathrm{d}\mathit{,}\mathrm{clust}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}\mathit{z}}{{\mathit{\chi }}^{\mathrm{2}}}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}{\mathit{a}}^{\mathrm{2}}\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathit{j̅}\mathrm{\left(}{\mathit{\nu }}^{\mathrm{\prime }}\mathit{,z}\mathrm{\right)}{\mathit{P}}_{\mathit{j}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,z}\mathrm{)}\mathit{,}& & \end{array}$(34)where ${\mathit{P}}_{\mathit{j}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ is the 3D power spectrum of the emissivities and is defined as follows: $\begin{array}{ccc}⟨\mathit{\delta j}\mathrm{\left(}{k}\mathit{,\nu }\mathrm{\right)}\mathit{\delta j}\mathrm{\left(}{{k}}^{\mathrm{\prime }}\mathit{,}{\mathit{\nu }}^{\mathrm{\prime }}\mathrm{\right)}⟩\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}\mathit{j̅}\mathrm{\right(}\mathit{\nu }\mathrm{\left)}\mathit{j̅}\mathrm{\right(}{\mathit{\nu }}^{\mathrm{\prime }}\mathrm{\left)}{\mathit{P}}_{\mathit{j}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\right(}{k}\mathrm{\left)}{\mathit{\delta }}^{\mathrm{3}}\mathrm{\right(}{k}\mathrm{-}{{k}}^{\mathrm{\prime }}\mathrm{)}\mathit{.}& & \end{array}$(35)In the context of CIB anisotropy modelling, the simplest version of the so-called “halo model”, which provides one view of the large-scale structure of the Universe as clumps of dark matter, consists of equating P_j with the galaxy power spectrum P_gg. This is equivalent to assuming that the CIB is sourced equally by all galaxies, so that the spatial variations in the emissivities trace the galaxy number density, $\begin{array}{ccc}\mathit{\delta j}\mathit{/}\mathit{j̅}\mathrm{=}\mathit{\delta }{\mathit{n}}_{\mathrm{gal}}\mathit{/}\mathit{n̅}\mathrm{gal}\mathit{.}& & \end{array}$(36)The dark matter haloes are populated through a halo occupation distribution (HOD) prescription. Ultimately, ${\mathit{P}}_{\mathrm{gg}}^{}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}$ is written as the sum of the contributions of galaxies within a single dark matter halo (“1h”) and galaxies belonging to two different haloes (“2h”): ${\mathit{P}}_{\mathrm{gg}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{1}\mathrm{h}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{2}\mathrm{h}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\hspace{0.17em}\mathit{.}$(37)On large scales P_2h reduces to a constant bias (squared) times the linear theory power spectrum, while the 1-halo contribution encapsulates the non-linear distribution of matter.

In this paper, we developed two approaches for the modelling of CIB anisotropies.

• The first one (Sect. 5.4) is very simple, and takesadvantage of the accurate measurement of CIB anisotropies withPlanck and IRIS at large angular scales. As analternative to the HOD model forP_gg we use a constant bias model in which ${\mathit{P}}_{\mathrm{gg}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{=}{\mathit{b}}_{\mathrm{eff}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{lin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}$, where b_eff is a redshift- and scale-independent bias and P_lin(k) is the linear theory, dark-matter power spectrum.

• The second one (Sect. 5.5) uses the anisotropies at all angular scales, and takes advantage of the frequency coverage of our measurements, to constrain a halo model with a luminosity–mass dependence. As a matter of fact, the model described above, which assumes that emissivity density traces galaxy number density (Eq. (36)), implies that all galaxies contribute equally to the emissivity density, irrespective of the masses of their host haloes. It assumes that all galaxies have the same luminosity, which is a crude assumption, as the luminosity and the clustering strength are closely related to the mass of the host halo.

Before becoming fully immersed in the details of our model building and fitting, we want to stress that the purpose of the following sections is to build as “physical” a model as possible that reproduces our measurements. To do so, we will expand upon the large amount of work in the literature that has exploited the simplifying concept of the halo model. We hope that in this way our work will have a wider impact. Nevertheless, it has to be emphasized that these models are very phenomenological and multi-faceted. To some extent, there is no such thing as the halo model, since there are many hidden assumptions underlying the application of the idea. For example, in our approach we will rely on a concentration prescription as a function of mass and redshift, and use this all the way to a redshift of 6, thus pushing into a regime where the model has not previously been validated. The same holds for our ansatz for the L–M relation and more generally to the concept of the HOD. To explore the dependence of our conclusions on these hidden assumptions is a task that goes well beyond the scope of the current paper and would certainly require a much more extensive use of simulations. It would of course be possible to weaken some of the assumptions by including them as Bayesian priors in our fit – but that would simply be side-stepping some of the serious conceptual limitations which exist here. For these reasons, we chose, in this current modelling effort, to simply make our assumptions clear and justify our choice of fixing values when possible. This approach is well-defined, but leaves a degree of uncertainty unaccounted for in the error budget, which should be kept in mind when interpreting the results, particularly in terms of physical parameters. Uncertainties in parameters would certainly be expected to increase if we relax some of the model’s assumptions.

5.3. Fitting for a model

The power spectra that are computed by the models need to be colour-corrected, from a CIB SED to our photometric convention νI_ν = const. We use the CIB SED from Béthermin et al. (2012a) to compute the colour corrections. They are equal to 1.076, 1.017, 1.119, 1.097, 1.068, 0.995, and 0.960 at 100, 143, 217, 353, 545, 857, and 3000 GHz, respectively. The correction to the power spectra follows ${\mathit{C}}_{\mathit{\ell ,\nu ,}{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{model}}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{cc}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{cc}}_{{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{=}{\mathit{C}}_{\mathit{\ell ,\nu ,}{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{measured}}\mathit{.}$(38)We use the same colour corrections (cc) for both the star-forming galaxy shot noise and the CIB power spectrum. For the radio galaxy shot noise we use a power law S_ν ∝ ν^α, with α = − 0.5. This is the average spectral index for radio sources that mainly contribute to the shot-noise power spectra. With such an SED, we find that the colour corrections are all lower than 0.7% for 100 ≤ ν ≤ 857 GHz. We thus neglect them.

To search for our best-fit model, we follow this scheme.

• 1.

  Take the residual map power spectra, as given in Table D.1.

• 2.

  Discard the first two bins at multipoles ℓ = 53 and 114 (following Sects. 3.2.2 and 3.2.3).

• 3.

  Correct for SZ-related residuals, ${\mathit{C}}_{\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$, and ${\mathit{C}}_{\mathrm{CIB}\mathrm{-}\mathrm{SZcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$, following Eqs. (8), (9), (12), and add the errors of these corrections quadratically to the error bars given in Table D.1.

• 4.

  Apply the colour correction to convert the theoretical model from measured Jy² sr^-1 to Jy² sr^-1 [νI_ν = const.] (Eq. (38)).

• 5.

  Compute the χ² value between the theoretical model and the observations, further applying the correction ${\mathit{C}}_{\mathrm{CIBcorr}}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}$ (Eq. (4)), and adding calibration errors, as described in the next item.

• 6.

  The calibration uncertainties are treated differently than the CIB power spectra error bars. We use an approach similar to the galaxy number counts model of Béthermin et al. (2011). A calibration factor f_cal is introduced. It has an initial value of 1, but can vary inside a Gaussian prior, centred on the calibration errors given in Table 1. We add a term to the χ² that takes into account the estimated calibration uncertainties: ${\mathit{\chi }}_{\mathrm{cal}}^{\mathrm{2}}\mathrm{=}{\sum }_{\mathrm{bands}}\mathrm{(}{\mathit{f}}_{\mathrm{cal}}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{2}}$, where ${\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{2}}$ are the calibration errors. The C_ℓ computed for the model is thus modified according to ${\mathit{C}}_{\mathit{\ell ,\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{c}\mathit{al}}\mathrm{=}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{cal}}{\mathit{f}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{cal}}{\mathit{C}}_{\mathit{\ell ,\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{model}}$.

5.4. Constraints on SFRD and effective bias from the large-angle linear scales

5.4.1. Fitting the linear model to the data

Planck is a unique probe of the large-scale anisotropies of the CIB. At ℓ ≲ 1000, the clustering is dominated by the correlation between dark matter haloes (the 2-halo term, e.g., Planck Collaboration 2011RPlanck Collaboration XVIII 2011). Planck data thus give the opportunity to put new constraints on both star-formation history and clustering of star-forming galaxies, using only the linear part of the power spectra. In this modelling, we consider only the 2-halo contribution to the cross-power spectrum between maps at frequency ν and ν′ (or auto-spectrum if ν = ν′), which can be written as ${\mathit{C}}_{\mathit{\ell ,\nu ,}{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{2}\mathrm{h}}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}\mathit{z}}{{\mathit{\chi }}^{\mathrm{2}}}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}{\mathit{a}}^{\mathrm{2}}{\mathit{b}}_{\mathrm{eff}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathit{j̅}\mathrm{\left(}{\mathit{\nu }}^{\mathrm{\prime }}\mathit{,z}\mathrm{\right)}{\mathit{P}}_{\mathrm{lin}}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,z}\mathrm{)}\mathit{,}$(39)and we fit only for ℓ ≤ 600. Here b_eff is the effective bias of infrared galaxies at a given redshift, i.e., the mean bias of dark matter haloes hosting infrared galaxies at a given redshift weighted by their contribution to the emissivities. This term implicitly takes into account the fact that more massive haloes are more clustered. The link between this simple approach and the HOD approach of Sect. 5.5 is discussed in Appendix C. We compute P_lin(k) using CAMB⁷.

The emissivities are derived from the star formation density ρ_SFR following (see Appendix B) $\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{SFR}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{S}}_{\mathit{\nu ,}\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{K}}\mathit{,}$(40)where K is the Kennicutt (1998) constant (SFR/L_IR = 1.7 × 10^-10M_⊙ yr^-1 for a Salpeter IMF) and S_ν,eff(z) the mean effective SED of all infrared galaxies at a given redshift. They are deduced from the Béthermin et al. (2012a) model (see Appendix B). These SEDs are a mixture of secularly star-forming galaxies and starburst galaxies. The dust temperature here increases with redshift following the measurements of Magdis et al. (2012) (see Sect. 6.4 for a discussion about the choice of the SED library for the modelling and the impact on the results).

There are degeneracies between the evolution of the bias and of the emissivities. In order to break them, we put some priors on the following quantities.

• The local infrared luminosity density,ρ_SFR(z = 0) = (1.95 ± 0.3) × 10^-2M_⊙ yr^-1 (Vaccari et al. 2010), converted using the H₀ value measured by Planck).

• The local bias of infrared galaxies, b = 0.84 ± 0.11 (Saunders et al. 1992), converted using σ₈ measured by Planck.

• The mean level of the CIB deduced from galaxy number counts, $\mathrm{12.}{\mathrm{6}}_{-1.7}^{\mathrm{+}\mathrm{8.3}}\hspace{0.17em}\mathrm{nW}\hspace{0.17em}{\mathrm{m}}^{-2}\hspace{0.17em}{\mathrm{sr}}^{-1}$ at 3000 GHz from Berta et al. (2011), $\mathrm{6.}{\mathrm{5}}_{-1.6}^{\mathrm{+}\mathrm{1.7}}\hspace{0.17em}\mathrm{nW}\hspace{0.17em}{\mathrm{m}}^{-2}\hspace{0.17em}{\mathrm{sr}}^{-1}$ at 857 GHz, and $\mathrm{2.}{\mathrm{1}}_{-0.6}^{\mathrm{+}\mathrm{0.7}}\hspace{0.17em}\mathrm{nW}\hspace{0.17em}{\mathrm{m}}^{-2}\hspace{0.17em}{\mathrm{sr}}^{-1}$ at 545 GHz from Béthermin et al. (2012c), and finally >0.27 nW m^-2 sr^-1 from Zemcov et al. (2010) at 353 GHz. These values are colour-corrected from PACS, SPIRE and SCUBA to Planck and IRAS, using the Béthermin et al. (2012a) model.

In this simple analysis, we want to measure only two quantities: the effective bias and its evolution with redshift, b_eff(z); and the star formation density history, ρ_SFR(z).

Inspired by the redshift evolution of the dark matter halo bias, we chose the following simple parametric form for the evolution of the effective bias: ${\mathit{b}}_{\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{b}}_{\mathrm{0}}\mathrm{+}{\mathit{b}}_{\mathrm{1}}\mathit{z}\mathrm{+}{\mathit{b}}_{\mathrm{2}}{\mathit{z}}^{\mathrm{2}}\mathit{.}$(41)For the star formation history, the values of ρ_SFR at z = 0, 1, 2 and 4 are free parameters, and we connect these points assuming a power law in (1 + z), using the two last points to extrapolate ρ_SFR at z> 4. We perform a Monte Carlo Markov chain analysis of the global parameter space. We assume Gaussian uncorrelated error bars for uncertainties, which are a linear combination of statistical and beam errors. The calibration uncertainties are treated following the method described Sect. 5.2.

To be independent of the exact level of the Poisson and 1-halo power spectrum in our linear analysis, we fit only for ℓ ≤ 600 measurements. For such ℓs, contamination by the Poisson and 1-halo terms is lower than ~10% (except for 3000×3000 where it reaches ~25% at ℓ = 502, see Fig. 10). We nevertheless add to our model the small correction due to the 1-halo and Poisson terms, as derived from the Béthermin et al. (2013) model.

5.4.2. Results

With a best-fit χ² value of 35 for 41 degrees of freedom, we obtain a very good fit to the data. In Table 8 we quote median values and marginalized limits for the parameters. The posterior value of parameters for which we imposed a Gaussian prior (local effective bias and SFRD, plus calibration factors) are all within the 1σ range of the prior values (except the 857 GHz calibration factor which is at 1.2σ).

Figure 11 shows the evolution of the star formation density with redshift (upper panel). Our derived star formation history nicely agrees with the infrared measurements of the dust-obscured star-formation rate density of Rodighiero et al. (2010) and Magnelli et al. (2011), up to z ~ 2. At higher redshift, our determination is marginally compatible (2σ) with Gruppioni et al. (2013), but in very good agreement with the recent work of Burgarella et al. (2013) at z = 3. We also compared our measurements with the UV estimate of star formation (not corrected for dust-attenuation) from Bouwens et al. (2007), Cucciati et al. (2012), and Reddy & Steidel (2009). Below z ~ 3, the bulk of the UV light emitted by young, short-lived stars is reprocessed in the infrared. Above this redshift, we find that the star formation probed in the UV and IR regimes have roughly an equal contribution. The infrared regime alone is thus no longer a good measure of the total star-formation rate density.

We also studied the evolution of the effective bias (lower panel of Fig. 11). We measure an increase of the bias with redshift. In Fig. 11 we compare the evolution of the galaxy dark matter bias with that of dark matter haloes of various mass (from Tinker et al. 2008). Our results are compatible with the track of dark matter haloes with 1−3 10¹²M_⊙, corresponding to the halo mass of maximal efficiency of star formation, as found in recent works (e.g., Béthermin et al. 2012b; Wang et al. 2013; and Behroozi et al. 2013a; and compatible with the related lensing magnification study of Hildebrandt et al. 2013).

5.5. Halo model for CIB anisotropies

The halo model is a standard approach to describe the clustering of matter at all scales (Cooray & Sheth 2002). Starting from the assumption that all galaxies live in dark matter haloes, the clustering power spectrum can be considered as the sum of two components: the 1-halo term (labelled P_1h), due to correlations of galaxies within the same halo, is responsible for the small-scale clustering; while the 2-halo term (P_2h), sourced by galaxy correlations in different haloes, describes the large-scale clustering.

The galaxy power spectrum is completely characterized by four main ingredients: the halo bias between dark matter and haloes; the halo density profile, describing the spatial distribution of dark matter inside a given halo; the halo mass function, specifying the number density of haloes with a given mass; and a prescription for filling dark matter haloes with galaxies, the so-called Halo Occupation Distribution (HOD).

A common assumption in the simplest versions of the halo model is that all galaxies have the same luminosity, regardless of their host dark matter halo (Viero et al. 2009; Amblard et al. 2011; Planck Collaboration 2011RPlanck Collaboration XVIII 2011; Xia et al. 2012; Viero et al. 2013b). However, as has already been pointed out in Shang et al. (2012), both galaxy clustering and galaxy luminosity are linked to host halo mass so that, in a statistical way, galaxies situated in more massive haloes have more stellar mass and are more luminous. The lack of such a link between galaxy luminosity and host halo mass in the model can lead to an interpretation of the clustering signal on small angular scales being due to a significant overabundance of satellite haloes (as in Amblard et al. 2011) with respect to what is found in numerical simulations (see discussion in Shang et al. 2012; Viero et al. 2013b). However, this can instead be due to a smaller number of galaxies, but with higher luminosity.

In this paper, we assume a halo model with a galaxy luminosity-halo mass relation similar to the one introduced in Shang et al. (2012) and also used in Viero et al. (2013b). We define haloes as overdense regions with a mean density equal to 200 times the mean density of the Universe and we assume an NFW profile (Navarro et al. 1997) for the halo density profile, with a concentration parameter as in Cooray & Sheth (2002). Fitting functions of Tinker et al. (2008) and the associated prescription for the halo bias (see Tinker et al. 2010) will be used for the halo and subhalo mass functions, respectively.

In the next subsections we will introduce the halo model that we use and we will describe how our analysis constrains its main parameters.

5.5.1. A halo model with luminosity dependence

The relation between the observed flux S_ν and the luminosity of a source at a comoving distance χ(z) is given by: $\begin{array}{ccc}{\mathit{S}}_{\mathit{\nu }}& \mathrm{=}& \frac{{\mathit{L}}_{\mathit{\nu }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathrm{4}\mathit{\pi }{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathit{,}\end{array}$(42)and the galaxy emissivity can be written as $\begin{array}{ccc}\mathit{j̅}\mathit{\nu }\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \end{array}$(43)where dn/dL denotes the infrared galaxy luminosity function.

In general, in order to model the galaxy luminosity-halo mass relation, we should introduce a scatter describing the probability density P(L | M) for a halo (or a subhalo) of mass M to host a galaxy with luminosity L (as in the conditional luminosity function models of, e.g., Yang et al. 2003, 2005; Cooray & Milosavljević 2005; Cooray 2006; Amblard & Cooray 2007; and De Bernardis & Cooray 2012). In order to keep the analysis as simple as possible, we neglect any scatter and introduce L_{cen,ν(1 + z)}(M_H,z) (for central galaxies) and L_{sat,ν(1 + z)}(m_SH,z) (for satellite galaxies), where M_H and m_SH denote the halo and subhalo masses, respectively, Eq. (43) can be re-written as: $\begin{array}{ccc}\mathit{j̅}\mathit{\nu }\mathrm{\left(}\mathit{z}\mathrm{\right)}& \mathrm{=}& \int \mathrm{d}\mathit{M}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\{\frac{}{}{\mathit{N}}_{\mathrm{cen}}{\mathit{L}}_{\mathrm{cen}\mathit{,}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{\nu }}\mathrm{\left(}{\mathit{M}}_{\mathrm{H}}\mathit{,z}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{\int }\mathrm{d}{\mathit{m}}_{\mathrm{SH}}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{m}}\mathrm{\left(}{\mathit{m}}_{\mathrm{SH}}\mathit{,z}\mathrm{\right)}{\mathit{L}}_{\mathrm{sat}\mathit{,}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{\nu }}\mathrm{\left(}{\mathit{m}}_{\mathrm{SH}}\mathit{,z}\mathrm{\right)}\}\mathit{,}\end{array}$(44)where dn/ dm denotes the subhalo mass function and N_cen is the number of central galaxies inside a halo.

Introducing ${\mathit{f}}_{\mathit{\nu }}^{\mathrm{cen}}$ and ${\mathit{f}}_{\mathit{\nu }}^{\mathrm{sat}}$ for central and satellite galaxies, $\begin{array}{ccc}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{cen}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}{\mathit{N}}_{\mathrm{cen}}\frac{{\mathit{L}}_{\mathrm{cen}\mathit{,}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{\nu }}\mathrm{\left(}{\mathit{M}}_{\mathrm{H}}\mathit{,z}\mathrm{\right)}}{\mathrm{4}\mathit{\pi }}\mathit{,}& & \end{array}$(45)$\begin{array}{ccc}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{sat}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathit{M}\\ {\mathit{M}}_{\min }\end{array}\mathrm{d}\mathit{m}\frac{\mathrm{d}{\mathit{n}}_{}}{\mathrm{d}\mathit{m}}\mathrm{\left(}{\mathit{m}}_{\mathrm{SH}}\mathit{,z}\mathrm{\right|}\mathit{M}\mathrm{)}\\ & & \end{array}$(46)then the power spectrum of CIB anisotropies at observed frequencies ν,ν′ can be written as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{1}\mathrm{h}\mathit{,\nu }{\mathit{\nu }}^{\mathrm{\prime }}}^{}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{j̅}\mathit{\nu }\mathit{j̅}{\mathit{\nu }}^{\mathrm{\prime }}}{\mathrm{\int }}_{{\mathit{M}}_{\min }}^{\mathrm{\infty }}\mathrm{d}\mathit{M}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{M}}\\ & & \mathrm{\times }\hspace{0.17em}\mathrm{\{}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{cen}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}{\mathit{f}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{sat}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathit{u}\mathrm{(}\mathit{k,M,z}{\mathrm{)}}^{}\\ & & \mathrm{+}{\mathit{f}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{cen}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{sat}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathit{u}\mathrm{\left(}\mathit{k,M,z}\mathrm{\right)}\\ & & \mathrm{+}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{sat}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}{\mathit{f}}_{{\mathit{\nu }}^{\mathrm{\prime }}}^{\mathrm{sat}}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}{\mathit{u}}^{\mathrm{2}}\mathrm{(}\mathit{k,M,z}{\mathrm{)}}^{\mathrm{\}}}\mathit{,}\\ {\mathit{P}}_{\mathrm{2}\mathrm{h}\mathit{,\nu }{\mathit{\nu }}^{\mathrm{\prime }}}^{}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{j̅}\mathit{\nu }\mathit{j̅}{\mathit{\nu }}^{\mathrm{\prime }}}{\mathit{D}}_{\mathit{\nu }}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}{\mathit{D}}_{{\mathit{\nu }}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}{\mathit{P}}_{\mathrm{lin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathit{.}\end{array}$Here $\begin{array}{ccc}{\mathit{D}}_{\mathit{\nu }}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{c}\mathrm{\infty }\\ {\mathit{M}}_{\min }\end{array}\mathrm{d}\mathit{M}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{M}}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathit{u}\mathrm{\left(}\mathit{k,M,z}\mathrm{\right)}\\ & & \end{array}$(49)with u(k,M,z) being the Fourier transform of the halo density profile.

5.5.2. Parameterizing the L–M relation

In the simplest version of the halo model, where galaxies residing in haloes of different masses have the same luminosity, the galaxy power spectrum is fully determined by the HOD, namely the function describing the number of central and satellite galaxies in each dark matter halo. In the model used here, the power spectrum depends, additionally, on the function L_{(1 + z)ν}(M_H,z), where M denotes the halo mass. The luminosity L_{(1 + z)ν}(M_H,z) depends on three variables: the redshift z; the mass of the host (sub)halo; and the observing frequency ν. We will consider the following assumptions about the structure of the luminosity–mass relation L(M).

• We assume no difference between haloes and subhaloes with the same mass, so that L(M_H,z) = L(m_SH,z), for M_H = m_SH. While recent studies (e.g., Rodríguez-Puebla et al. 2012, 2013) show some indication that satellite galaxies tend to have slightly more stellar mass than central galaxies with the same halo mass, these results depend on the subhalo mass definition used; in particular, the luminosity–mass relation for satellites and central galaxies has been found to be not very different when the mass of the subhalo is defined at the time of accretion (as done in this paper).

• A very simple functional form (see Blain et al. 2003, and reference therein) is assumed for galaxy SEDs: $\begin{array}{ccc}\mathrm{\Theta }\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{\propto }\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$(50)Here B_ν denotes the Planck function, while the emissivity index β gives information about the physical nature of dust and in general depends on grain composition, temperature distribution of tunnelling states and wavelength-dependent excitation (e.g., Meny et al. 2007). The power-law function is used to temper the exponential Wien tail at high frequencies and obtain a shallower SED shape, more in agreement with observations. The temperature is assumed to be a function of redshift according to $\begin{array}{ccc}{\mathit{T}}_{\mathrm{d}}\mathrm{\equiv }{\mathit{T}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathit{\alpha }}\mathit{.}& & \end{array}$(51)This dependence of the temperature with redshift can be due to different physical processes, such as more compact geometries for galaxies at high redshift (Magdis et al. 2012), a global evolution of the SED (e.g., Addison et al. 2013; Béthermin et al. 2013) or the increase of the CMB temperature with redshift (Blain 1999). The SED functions at high and low frequencies are connected smoothly at the frequency ν₀ satisfying $\begin{array}{ccc}\frac{\mathrm{d}\ln \mathrm{\Theta }\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}}{\mathrm{d}\ln \mathit{\nu }}\mathrm{=}\mathrm{-}\mathit{\gamma }\mathrm{·}& & \end{array}$(52)The range of variation for the parameters α, γ, and ν_∗ is large enough to ensure that we do not exclude non-negligible regions of the multidimensional parameter space; however we assume physically motivated priors for both the temperature (T₀ in the range 20–60 K, see measurements in e.g., Dunne et al. 2000; Chapman et al. 2005; Amblard et al. 2010; and Hwang et al. 2010) and the emissivity index (β in the range 1.5–2.0). The correct choice of β is a matter of debate; measurements of Milky Way dust, and in external galaxies (e.g., Boselli et al. 2012), give values in the range 1–2, but allowing for some degree of correlation between dust temperature and emissivity index (see, e.g., Paradis et al. 2010) it is possible to obtain β> 2 for low dust temperatures (T_d ≤ 18 K). On the theoretical side, while models for both insulating and conducting materials naturally give β = 2 at long wavelengths (e.g., Draine & Lee 1984), significant deviations from the value β = 2 occur when accounting for the disordered structure of the amorphous dust grains (Meny et al. 2007). Indeed, some authors (Shang et al. 2012; Viero et al. 2013b) allow for values β> 2 when fitting CIB data. In this analysis we prefer to be conservative and, since we assume the condition T_d> 20 K, we also impose β ≤ 2; this will allow us to draw solid conclusions on the other parameters of the model, avoiding regions of the parameter space whose physical interpretation is questionable.

• We assume a redshift-dependent, global normalization of the L–M relation of the form $\begin{array}{ccc}\mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{z}}^{\mathrm{)}}}^{\mathit{\delta }}\mathit{.}& & \end{array}$(53)The parameter δ will be allowed to vary in the range 0–7. Such a redshift dependence can be justified considering the evolution of the specific far infrared luminosity (L_IR/M_∗) with redshift: if the ratio of stellar mass to halo mass evolves only mildly with redshift (see e.g., Neistein et al. 2011), then the ratio L_IR/M_H should evolve approximately as the specific infrared luminosity. The semi-analytic galaxy formation model of De Lucia & Blaizot (2007) shows the evolution of such a quantity with redshift as a power law with a slope of about 2.5, while observations performed by Oliver et al. (2010) indicate a much steeper slope, around 4.4.

• We assume a log-normal function for the dependence of the galaxy luminosity on halo mass: $\begin{array}{ccc}\mathrm{\Sigma }\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}\mathit{M}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathit{\sigma }}_{\mathit{L}\mathit{/}\mathit{M}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{\left(}{\log }_{\mathrm{10}}\mathrm{\right(}\mathit{M}\mathrm{)}\hspace{0.17em}\mathrm{-}\hspace{0.17em}{\log }_{\mathrm{10}}\mathrm{(}{\mathit{M}}_{\mathrm{eff}}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\sigma }}_{\mathit{L}\mathit{/}\mathit{M}}^{\mathrm{2}}}\mathit{.}& & \end{array}$(54)Here M_eff describes the peak of the specific IR emissivity, while the parameter σ_L/M describes the range of halo masses used for producing the IR luminosity; we will assume that ${\mathit{\sigma }}_{\mathit{L}\mathit{/}\mathit{M}}^{\mathrm{2}}\mathrm{=}\mathrm{0.5}$ throughout this paper and we checked that results do not significantly change when assuming ${\mathit{\sigma }}_{\mathit{L}\mathit{/}\mathit{M}}^{\mathrm{2}}\mathrm{=}\mathrm{0.65}$, as in Béthermin et al. (2012b). The reason for choosing a log-normal functional form is that star formation is active only over a given range of halo masses, being suppressed at both the low- and the high-mass end by mechanisms such as photoionization, supernovae heating, feedback from active galactic nuclei and virial shocks (see e.g., Benson et al. 2003; Croton et al. 2006); it is then possible to identify a peak in the L–M relation, which describes the maximum in the average infrared emissivity per unit mass.

• At the low-mass end, we assume a minimum mass M_min, which is a free parameter in the range 10¹⁰–10¹¹M_⊙, and we assume L = 0 for M<M_min.

The equation for the luminosity–mass relation can finally be written as $\begin{array}{ccc}{\mathit{L}}_{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{\nu }}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{=}{\mathit{L}}_{\mathrm{0}}\mathrm{\Phi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\Sigma }\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{\Theta }\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{\left)}\mathit{\nu }\mathrm{\right]}\mathit{,}& & \end{array}$(55)where L₀ is a free normalization parameter (which being not physically meaningful will not be further discussed).

5.5.3. Method and data used

In order to constrain the main parameters of our model, we fit for a total of 121 data points of the 15 possible combinations of Planck auto- and cross-power spectra at 217, 353, 545, 857, and 3000  GHz, considering the multipole range 187 ≤ ℓ ≤ 2649. We use the same procedure as described in Sect. 5.3 in order to colour-correct our model to the photometric convention νI_ν = const. and to include the corrections due to CIB over-subtraction and SZ-related residuals (see points 1–6 in Sect. 5.3). There are two main differences with respect to the linear model analysis outlined above (Sect. 5.4):

• we keep all the calibration parameters fixed atf^cal = 1, which assumption is justified by the analysis using the linear model, allowing us not to deal with too many parameters;

• we assume the same prior on the star formation rate density (Vaccari et al. 2010) as in the linear model but we do not use any constraints on the bias at redshift zero. We also assume flat priors on the mean level of the CIB at 545 and 857  GHz.

We perform a Monte Carlo Markov chain analysis of the global parameter space using a modification of the publicly available code CosmoMC (Lewis & Bridle 2002). We consider variations in the following set of eight halo model parameters: $\begin{array}{ccc}P\mathrm{\equiv }\mathrm{\left\{}\mathit{\alpha ,\beta ,\gamma ,\delta ,}{\mathit{M}}_{\mathrm{eff}}\mathit{,}{\mathit{M}}_{\min }\mathit{,}{\mathit{T}}_{\mathrm{0}}\mathit{,}{\mathit{L}}_{\mathrm{0}}\mathrm{\right\}}\mathit{.}& & \end{array}$(56)We assume the shot-noise levels given by the sum of the values quoted in Tables 6 and 7, from Béthermin et al. (2012a) and Tucci et al. (2011), respectively, and we assume flat priors around them with width given by their 1σ error.

The total number of free parameters in our analysis is then 23, consisting of the sum of eight halo model parameters plus 15 shot-noise parameters. While it is tempting to fix the shot-noise power spectra to their theoretically modelled values (in order not to deal with too many parameters and keep the analysis as simple as possible), we believe that, since these values are not very tightly constrained by their underlying models, it is better to let them vary as free parameters around their best estimates.

5.5.4. Results

With a best-fit χ² of 100.7 and 98 degrees of freedom, we obtain a remarkably good fit to the data. In Table 9 we quote mean values and marginalized limits for the model parameters. In the following, we comment on the results obtained for some parameters of the model and for some derived quantities.

Peak mass M_eff –

The mean value of the most efficient halo mass for generating the CIB, log (M_eff/M_⊙) = 12.6 ± 0.1, is in good agreement with results obtained from a similar analysis using Herschel CIB data at 250, 350 and 500 μm (Viero et al. 2013b), and with other analyses, using previous Planck and Herschel data (e.g., Shang et al. 2012; Xia et al. 2012), while it is slightly higher than results from other observations and simulations (e.g., Moster et al. 2010; Behroozi et al. 2013b; Béthermin et al. 2012b; Wang et al. 2013). We also checked for a possible redshift evolution of M_eff (which can be justified in the framework of the so-called “downsizing” idea), performing an MCMC run with the functional form ${\mathit{M}}_{\mathrm{eff}}\mathrm{=}{\mathit{M}}_{\mathrm{0}}\mathrm{\times }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathit{q}}\mathit{.}$(57)The large degeneracy between M₀ and q leads to very high values of M₀. The bias and SFRD have the same redshift evolution as in the case q = 0, but with much larger error bars (they are a factor of 6 times higher for the bias, for example).

Constraints on the dust temperature –

Parameterizing the average dust temperature of sources as ${\mathit{T}}_{\mathrm{d}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathit{\alpha }}\mathit{,}$(58)the data suggest a redshift evolution of the temperature, with T₀ = (24.4 ± 1.9) K and α = 0.36 ± 0.05. Such a trend, implying some kind of SED evolution, has been also found in e.g., Addison et al. (2013); Viero et al. (2013b). Experimental results from different surveys appear to have been quite contradictory, with systematics playing a critical role (e.g., Chapman et al. 2005; Coppin et al. 2008; Pascale et al. 2009; Amblard et al. 2010; Hwang et al. 2010; Chapin et al. 2011). But recently, some consensus has emerged on a scenario with an increase of dust temperature with redshift (Magdis et al. 2012; Viero et al. 2013a). The increase of temperature may be explained by a harder interstellar radiation field at earlier times (see Magdis et al. 2012, for a detailed discussion).

Constraints on the bias –

Galaxies are considered as a biased tracer of the dark matter field. The galaxy overdensity δ_g(k,z) is assumed to trace the underlying dark matter field δ_dm(k,z) via ${\mathit{\delta }}_{\mathrm{g}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{=}\mathit{b}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}{\mathit{\delta }}_{\mathrm{dm}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathit{,}$(59)where b(k,z) is the galaxy bias, which in general can depend not only on scale and redshift but also on luminosity, spectral type and colour. On large scales, the bias is generally assumed to be scale-independent; however, both numerical simulations (Kauffmann et al. 1999) and recent results from galaxy-galaxy lensing and galaxy clustering also indicate an increase of the bias with redshift (e.g., Mandelbaum et al. 2013), while Tegmark & Peebles (1998) show that the bias must be close to unity when approaching z = 0. The combination of CMB and large-scale clustering data yields a bias parameter b ~ 1 (Verde et al. 2002) while Saunders et al. (1992) found bσ₈ = 0.84 ± 0.11 for IRAS galaxies, which, assuming σ₈ = 0.8, gives b ~ 0.86. In Fig. 11 we show our estimate of the redshift dependent bias; it is remarkable that, without assuming any prior on the value of the bias at redshift zero, we are able to obtain a very good fit to observations, with b(z = 0) = 1.1 ± 0.02.

Constraints on star formation history –

The mean value and 68% CL bounds on the cosmic star formation rate density ρ_SFR are plotted in Fig. 11. The parameter ρ_SFR has been computed following Eq. (40), replacing s_ν,eff(z) by the halo model SED as given in Eq. (50). We use the Kennicutt (1998) constant to convert infrared luminosity to star formation rate (SFR/L_IR = 1.7 × 10^-10M_⊙ yr^-1 for a Salpeter IMF). As can be seen from Fig. 11, the star formation rate densities predicted by both models used in this paper are in very good agreement for redshifts z ≲ 2 while there is a significant difference at higher redshifts. Interestingly, Behroozi et al. (2013b) reports a compilation of results, showing how measurements performed before 2006 predict quite a high value for the high redshift SFRD (then more compatible with the halo model results, see also the compilation in Hopkins & Beacom 2006), while measurements performed after 2006 and obtained with different assumptions about the dust present at z> 3, show a rapid decrease of the SFRD, which then becomes more compatible with results obtained using the linear model. In an attempt to reproduce the break in the SFRD seen in Fig. 14 for the linear model, we also imposed the condition δ = 0 at redshift z = 2 in the redshift normalization parameter Φ(z) of the L–M relation. Such a condition has already been considered in Shang et al. (2012) and, while it is motivated by some observations, it is also hard to explain from a theoretical point of view (e.g., Weinmann et al. 2011). In this case, we are able to obtain lower values of the SFRD at high redshifts (see Fig. 13), but at the price of slightly degrading the quality of the fit. We finally note that another potential reason for the discrepancy found at high redshift can be the difference in the inferred bias evolution between the linear model and the halo model: in fact, a lower value of the bias at high redshift (as found in the context of the halo model) can be compensated by higher values for the SFRD in the same redshift range. We will come back to the SFRD discussion in Sect. 6.4.

Finally, we are also able to determine the mean CIB intensity at the Planck frequencies considered in the analysis. The values obtained are presented in Table 10. Although higher, they are compatible, within the 95% CL, with results obtained from number counts measurements (Béthermin et al. 2012c).

5.6. CIB-CMB lensing cross-correlation

We tested the validity of our approach by comparing the predictions for our best-fit models with the measurements of the cross-correlation between the CIB and the CMB lensing potential presented in Planck Collaboration 2013RPlanck Collaboration XVIII (2014). For the linear model, we computed the cross-correlation following: ${\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu \phi }}\mathrm{=}\mathrm{\int }{\mathit{b}}_{\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{j̅}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\frac{\mathrm{3}}{{\mathit{\ell }}^{\mathrm{2}}}{\mathrm{\Omega }}_{\mathit{rmm}}\mathit{,}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}\left(\frac{{\mathit{\chi }}_{\mathrm{\ast }}\mathrm{-}\mathit{\chi }}{{\mathit{\chi }}_{\mathrm{\ast }}\mathit{\chi }}\right){\mathit{P}}_{\mathrm{lin}}\mathrm{(}\mathit{k}\mathrm{=}\mathit{\ell }\mathit{/}\mathit{\chi ,z}\mathrm{)}\mathrm{d}\mathit{\chi ,}$(60)where χ_∗ is the comoving distance to the CMB last-scattering surface. We use a similar equation for the extended halo model.

Figure 15 shows a comparison between the model and the data. The halo model (as well as its variant with a break in the temperature and global normalization of the L–M relation at redshift z ~ 4, see Sect. 6.4) agrees remarkably well with the measurements for all channels. The linear model gives a higher prediction at 217 and 353  GHz (although compatible with the data points at the 1σ level).

6. Discussion

6.1. The 143 GHz case

Removing the CMB anisotropy at 143  GHz is very problematic, since the CMB power spectrum is about 5000 higher than the CIB at ℓ = 100. However, thanks to the exceptional quality of the Planck data, and the accuracy of the 100 and 143  GHz relative photometric calibration, we can obtain interesting measurements, using the same method to clean the maps and measure the power spectra as for the other channels. We show the measurements in Fig. 16, together with the best-fit CIB model. This estimate has been obtained by correcting the measurements for the SZ and spurious CIB (induced by the use of the 100  GHz map as a CMB template). Those corrections are important, especially for power spectra at low frequencies. For example, for 143 × 217, the SZ-related corrections decrease the measurements by 10−20%, while the correction for the spurious CIB (including the shot noise) increases the measurements by 30−60%. Since these corrections are large, we have not attempted to include the 143  GHz measurements when constraining the model, and we recommend the reader to take them with caution.

We show in Fig. 16 a comparison between extrapolation of the halo best-fit model to the 143  GHz cross-power spectra and our CIB power-spectrum estimates. The 143 × ν cross-power spectra agree quite well for ℓ < 1000, at least for ν ≥ 353 GHz. The 143 × 217 CIB power spectrum lies about 2σ above the prediction at intermediate scales (ℓ = 502 and 684). This CIB overestimate increases for the 143 × 143 power spectrum, which is certainly the most difficult to obtain; this is in excess with respect to the prediction for 300 <ℓ < 1000. At this frequency, however, the CIB auto-power spectrum measurements have to be taken with caution, as the correction for the spurious CIB can be as high as 70%, and is thus highly model dependent.

6.2. Frequency decorrelation

Using the power spectrum measurements, we can quantify the frequency decoherence. We measure the correlation between bands by averaging the quantity ${\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}\mathit{/}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathit{\nu }\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathit{\nu }}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{C}}_{\mathit{\ell }}^{{\mathit{\nu }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathit{\nu }}^{\mathrm{\prime }}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ for 150 <ℓ < 1000. We restrict ourselves to this ℓ range to have only the clustered CIB contribution (not the shot noise). Results are given in Table 11. We see that the CIB for the four HFI frequencies (from 217 to 857  GHz) is very well correlated, the worst case being between the 857 and 217  GHz channels, with a correlation of about 0.85. On the other hand, the correlation of all HFI bands with IRIS is quite low, between about 0.2 and 0.32 (with a large dispersion). This is expected, because the redshift distribution of CIB anisotropies evolves strongly between 3000 and ≤857 GHz, being biased towards higher redshifts at lower frequencies (e.g., Béthermin et al. 2013).

Contrary to the range 217–857  GHz the band correlation strongly varies with ℓ at 143  GHz, decreasing from ℓ = 150 to 1000 (see Table 11). Such a decrease is expected, based on the high shot-noise contribution at this frequency (see Fig. 16), and the fact that the shot noise is dominated by the radio source contribution, contrary to the ≥217  GHz channels. Moreover, for dusty star-forming galaxies, Béthermin et al. (2013) observed that the band correlation is lower for the shot noise than for correlated anisotropies; this also mimics a scale dependence when the shot-noise contribution is important.

6.3. Comparison with recent measurements

We now compare the CIB auto-spectrum measurements with the most recent measurements from Herschel-SPIRE (Viero et al. 2013b) and the earlier measurements from Planck (Planck Collaboration 2011RPlanck Collaboration XVIII 2011). We compute the SPIRE-HFI colour corrections using the CIB SED from Gispert et al. (2000) and the most recent bandpasses (see Planck Collaboration 2013IPlanck Collaboration IX 2014). In order to compare with HFI at 857 and 545  GHz, the power spectra at 350 and 500 μm have to be multiplied by 1.016 and 0.805, respectively. We use the SPIRE power spectra with only extended sources masked (such that the Poisson contribution is the same in both measurements). We see from Fig. 17 that the agreement between the SPIRE and HFI measurements (red circles versus blue circles) is excellent at 857  GHz. At 545  GHz, although compatible within the error bars there is a small difference, with the SPIRE power spectrum being higher than HFI by about 7% for 650 <ℓ < 1800 and by about 30% for 200 <ℓ < 600.

Between the publication of Planck Collaboration 2011RPlanck Collaboration XVIII (2011) and this paper, the photometric calibration of the two high-frequency HFI channels has been modified (see Planck Collaboration 2013HPlanck Collaboration VIII 2014). Using a planet-based calibration rather than a FIRAS-based calibration leads to a division of the calibration factors by 1.07 and 1.15 at 857 and 545  GHz, respectively. After correcting for these factors, the two Planck CIB measurements agree within 1σ at 545  GHz  and within 2σ at 857  GHz. At 217  GHz, the discrepancy we observe between the earlier Planck CIB measurements and those presented in this paper is explained by the SZ and CIB contamination of the 143  GHz-based CMB template used in Planck Collaboration 2011RPlanck Collaboration XVIII (2011). Note that with the new measurements, we do not improve the error bars, since the use of the 100  GHz channel as a CMB template adds more noise than the use of the 143  GHz channel. The apparent difference in shape between the CIB at 217  GHz and at higher frequencies has to be attributed to the shot noise, whose contribution relative to the correlated part is higher at 217  GHz, making the measured CIB flatter at this frequency.

For cross-power spectra, we can compare our determination with Hajian et al. (2012). We show in Fig. 18 the 857 × 217 and 857 × 143Planck power spectra, and those obtained from the cross-correlation between BLAST and ACT data. For this comparison, the shot noise contributions have been removed, as they are very different for BLAST, ACT and Planck. We have not applied any colour correction. Even if it appears that three points overlap in scale, the comparison can only be done at one of them, since the low-ℓ measurements from BLAST×ACT are very noisy. For this ℓ = 1750 scale, the two measurements agree within 1σ. This plot illustrates the complementarity between Planck and the high-ℓ measurements from ACT and SPT – interesting constraints will come from combining the various ℓ ranges and frequencies probed. However, a full analysis will require careful consideration of the different instrumental bandpasses, as well as the treatment of foregrounds and point source masking. At a first look, we can compare the CIB band-powers in the three experiments by extrapolating our best-fit model up to ℓ = 3000. We obtain and $\mathrm{5.9}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}_{\mathrm{CMB}}^{\mathrm{2}}$ at 217 and 143  GHz, respectively. These numbers have to be compared with and $\mathrm{(}\mathrm{4.8}\mathrm{\pm }\mathrm{0.9}\mathrm{)}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}_{\mathrm{CMB}}^{\mathrm{2}}$ for ACT at 219.6 and 149.7  GHz, respectively (Dunkley et al. 2013); and and $\mathrm{(}\mathrm{6.2}\mathrm{\pm }\mathrm{0.5}\mathrm{)}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}_{\mathrm{CMB}}^{\mathrm{2}}$ for SPT at 219.6 and 153.8  GHz, respectively (Reichardt et al. 2012).

6.4. The history of star formation density

The star formation histories recovered from the two different modelling approaches presented in Sects. 5.4 and 5.5 are consistent below z = 2, and agree with recent estimates of the obscured star-formation density measured by Spitzer and Herschel. At higher redshift, there are discrepancies between our two models and the estimate of Gruppioni et al. (2013). The linear model is about 1 and 2σ lower than their measurements at z = 2.5 and 3.5, respectively, while the halo model lies about >3σ above these data points. Such estimates assume a shape for the infrared luminosity function. They are strongly dependent on the faint-end slope assumption, since no data are available below the break of the luminosity function. This shows how measurement of the obscured star formation rate density at z> 3 is difficult.

We investigated the origin of the discrepancy between the two modelling approaches. In particular, we modified the halo model in two ways to see how different assumptions on the parametrization of the model can affect the results.

• 1.

  We fit the data by imposing the condition δ = 0 for redshifts z ≥ z_break = 2 in the redshift normalization parameter Φ(z) of the luminosity–mass relation. This parametrization, although degrading the quality of the fit somewhat, decreases the SFRD by a factor of about 5 for z = 4, which is now compatible with the linear model SFRD (see Fig. 13).

• 2.

  We fit the data by imposing two conditions: δ = 0 for redshift z ≥ z_break and T(z ≥ z_break) = constant = T(z_break) for z_break in the range 2–5. We find z_break = 4.2 ± 0.5, as can be seen in Fig. 13. In this case the SFRD is only reduced at very high redshift, by a factor 3 at z = 6. Note that allowing for a redshift break in the redshift evolution of the temperature avoids reaching unphysically high values at very high redshift.

We also show in Fig. 13 the SFRD measurements from the CIB-CMB lensing cross-correlation (Planck Collaboration 2013RPlanck Collaboration XVIII 2014). This compares favourably with a high SFRD level at high redshift.

Part of the discrepancy between the two modelling approaches can also be attributed to the effective bias. A higher bias, as that recovered at high redshift from the linear model, favours a lower SFRD.

We finally compared the SEDs used in the two approaches. The effective SEDs present a broader peak than the extended halo model SEDs, because they take into account the dispersion in dust temperature and the mixing between secularly star-forming galaxies and episodic starbursts. At z> 2, there are discrepancies in the Rayleigh-Jeans regime between the two templates, the effective SEDs being higher than the extended halo model SEDs; this discrepancy increases with redshift. At z = 5, for the same L_IR, the parametric SEDs of the halo model emits about 3 times less infrared light than the effective SEDs (and hence about an order of magnitude less fluctuations). This explains why this model requires a much higher star formation rate density to fit CIB anisotropies than the linear model⁸. Following the parametrization of Eq. (50), we fit the effective SEDs with a modified black body with a ν^1.75 emissivity law to obtain the dust temperature. In Fig. 19, we show the redshift evolution of the temperature of the two templates. Compared to the recent average temperatures of star-forming galaxies found by Viero et al. (2013a) up to z ~ 4, the temperature of the effective SEDs is a bit low while the temperature of the SEDs of the extended halo model is a bit high. Knowing the SEDs of the galaxies that are responsible for the bulk of the CIB is the principal limitation in our modelling framework. Accurate future measurement of the SEDs will be crucial to properly estimate the obscured star formation rate density at high redshift from the CIB anisotropies. This is important if one wants to determine whether or not the bulk of the star formation is obscured at high redshift, and whether the UV and Lyman-break galaxy populations are a complete tracer of the star formation in the early Universe.

6.5. CIB non-Gaussianity

Lacasa et al. (2012) proposed a phenomenological prescription for the CIB bispectrum based on its power spectrum, namely ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\mathit{\alpha }\sqrt{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}\hspace{0.17em}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}\hspace{0.17em}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}}\mathit{,}$(61)where α is a dimensionless parameter quantifying the intrinsic level of non-Gaussianity. Using the best-fit power spectrum of the CIB model described in Sect. 5.5, we fitted this parameter α through a χ² minimization using the covariance matrix described in Sect. 4.3. The resulting best-fit α, its error bar (computed using a Fisher matrix analysis), and the χ² value of the best fit can be found in Table 12.

The consistency of α across frequencies shows that the measured bispectrum has a frequency dependence consistent with that of the power spectrum. The best-fit α values are consistent with the values predicted using the number counts model of Béthermin et al. (2011). The best-fit αs are of the same order, although a little lower than those found by Lacasa et al. (2012) on simulations by Sehgal et al. (2010), since they found α ≃ 3 × 10^-3. This indicates a lower level of CIB non-Gaussianity than in the simulations by Sehgal et al. (2010).

The χ² value of the fit shows that the prescription does not provide a very good model of the data as frequency increases; visual inspection reveals that this mainly comes from the fact that the measured bispectrum has a steeper slope than the prescription. To quantify the slope of the measured bispectrum, we fit a power law to the measurements, i.e., ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\mathit{A}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\left(\frac{{\mathit{\ell }}_{\mathrm{1}}\hspace{0.17em}{\mathit{\ell }}_{\mathrm{2}}\hspace{0.17em}{\mathit{\ell }}_{\mathrm{3}}}{{\mathit{\ell }}_{\mathrm{0}}^{\mathrm{3}}}\right)}^{\mathrm{-}\mathit{n}}\mathit{,}$(62)where we chose as the pivot scale ℓ₀ = 320, which is the centre of the second multipole bin. Table 13 presents the obtained best-fit values for the amplitude A and the index n, as well as their error bars and correlation (computed again with Fisher matrices) and the χ² value.

The power law provides a significantly better fit to the data than the prescription of Eq. (61), having lower best-fit χ² values. The indices obtained are coherent between frequencies, and significantly steeper than the Eq. (61) prescription, which is n ~ 0.6 (since C_ℓ ∝ ℓ^-1.2).

There is no sign of flattening of the bispectrum, showing that the shot-noise contribution is subdominant in this multipole range. This is consistent with shot-noise estimates based on the number counts model of Béthermin et al. (2011).

A detection of CIB non-Gaussianity has recently been reported by Crawford et al. (2014). In that paper, they used the prescription proposed by Lacasa et al. (2012) to give the amplitude of the bispectrum. In comparison with this analysis, we provide a higher detection significance, and at several frequencies. Most importantly, we find an indication that the CIB bispectrum is steeper than the prescription of Lacasa et al. (2012), although well fitted by a power law. However, our steeper CIB bispectrum at 217  GHz, extrapolated up to ℓ = 2000, is compatible with the Crawford et al. (2014) bispectrum measurement at ℓ = 2000 within 1σ. The steeper slope may be an indication that the contribution of more massive haloes to the CIB bispectrum is smaller than in the models studied by Lacasa et al. (2014) and Pénin et al. (2014).

7. Conclusions

We have presented new measurements of the CIB anisotropies with Planck. 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 have been able to measure the clustering of dusty, star-forming galaxies at 143, 217, 353, 545, and 857 GHz, with unprecedented precision. For the fist time we also measured the bispectrum from ℓ ≃ 130 to 900 at 217, 353, and 545 GHz. The CIB power spectrum is also measured with IRAS at 3000  GHz.

We worked on 11 independent fields, chosen to have high angular-resolution Hi data and low foreground contamination. The total areas used to compute the angular power spectrum is about 2240 deg². This improves over previous Planck and Herschel analyses by more than an order of magnitude. For the bispectrum, the total area is about 4400 deg².

To obtain the CIB, the HFI and IRAS maps were cleaned using two templates: Hi for Galactic cirrus; and the Planck 100  GHz map for CMB. We used new Hi data that covers very large portions of the sky. The large areas forced us to build a dust model that takes into account the submillimetre-Hi emissivity variations. However, because the Hi is not a perfect tracer of dust emission (e.g., the dark gas), and clearly contains dust-deficient clouds, we had to reduce the sky fraction to the lowest Hi column-density parts of the sky. The 100 GHz Planck channel, cleaned of Galactic dust and sources, and then filtered, provides a good template for the CMB. This is because it has an angular resolution close to the higher frequency channels, from which we measure the CIB, and has the advantage of being an “internal” template, meaning that its noise, data reduction processing steps, photometric calibration, and beam are all well known. It has the drawback of contaminating CIB measurements with tSZ signal and spurious CIB coming from the correlation between the CIB at 100 GHz, and the CIB at higher frequency. The tSZ and spurious CIB corrections are relatively small for frequencies ν ≥ 217 GHz. At 143  GHz, while the tSZ and tSZ–CIB corrections are still rather small (lower than 20%), the spurious CIB is a very large correction (between 30 and 70% at intermediate scales) due to the high level of CIB correlation between 100 and 143  GHz. Thus, the 143  GHz CIB measurements strongly rely on the CIB model used to compute the correction.

Due to dust contamination at high frequency, as well as radio sources, SZ and CIB at low frequency, we conservatively restrict our bispectrum measurement to the three frequencies 217, 353, and 545  GHz. We measure the bispectrum due to the clustering of dusty star forming galaxies from ℓ ≃ 130 to 900. It is detected with a very high significance, >28σ for all configurations, and even >4.5σ for individual configurations at 545  GHz. Such measurements are completely new; they open a window for constraining models of CIB source emission that we have not yet fully explored in this paper.

We developed two approaches for modelling the CIB anisotropies. The first takes advantage of the accurate measurement of CIB anisotropies performed with Planck and IRIS at large angular scales, and uses only the linear part of the power spectra. The second approach uses the measurements at all angular scales, and takes advantage of the frequency coverage, to constrain a halo model with a luminosity–mass dependence. We find that both models give a very good fit to the data. Our main findings are as follows.

• The models give strong constraints on the star formation history up to redshift ~2.5. At higher redshift, the accuracy of the star formation history measurement is strongly degraded by the uncertainty on the SED of CIB galaxies. An accurate measurement of SEDs of galaxies that are responsible for the bulk of the CIB will be crucial to estimate properly the obscured star-formation rate density at high redshifts from the CIB anisotropies.

• As found in other recent studies, haloes of mass M_eff = 4 × 10¹²M_⊙ appear to be the most efficient at actively forming stars.

• CIB galaxies have warmer temperatures as redshift increases (T_d(z) = 24.4 × (1 + z)^0.36 K in the extended halo model). This is compatible with the most recent Herschel observations, and can be explained by a harder interstellar radiation field in high-z galaxies.

• The same halo occupation distribution can simultaneously fit all power spectra. However, the 1-halo term is significantly reduced compared to previous studies (Pénin et al. 2012a; Planck Collaboration 2011RPlanck Collaboration XVIII 2011). This is due to a lower contribution to the clustering from low-z massive haloes, as also observed in Béthermin et al. (2013).

We find that the CIB bispectrum is steeper than the prescription developed by Lacasa et al. (2012). Just like the reduction of the 1-halo term in the power spectrum, this may be an indication that the contribution of massive haloes to the CIB bispectrum is smaller than in the models studied by Lacasa et al. (2014) and Pénin et al. (2014). The bispectrum is quite well fitted by a power law. This can be used to provide valuable constraints on the potential contamination of measurements of the primordial CMB bispectrum on large scales.

While our component separation process is successful in extracting the CIB from the maps, the next step is to use a full multi-frequency fitting procedure to separate the CIB power spectrum and bispectrum, from the tSZ (and kSZ) effects, the CMB, and the extragalactic source contribution. Simultaneously taking into account all the components will improve our ability to separate them. The goal is to give unprecedented limits on the reionization history of the Universe, as well as understanding the history of star formation in dark matter haloes.

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP-PNCG (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN and JA (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members with the technical or scientific activities they have been involved into, can be found at http://www.rssd.esa.int/index.php?project=PLANCK&page=PlanckCollaboration. The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. Some Hi data used in this paper are based on observations with the 100 m telescope of the MPIfR (Max-Planck-Institut für Radioastronomie) at Effelsberg.

References

Appendix A: H I data

We describe in this Appendix the three Hi surveys that we use to remove Galactic dust contamination from the frequency maps. The 21-cm Hi spectra were obtained with: (1) the 100-m Green Bank Telescope (GBT); (2) the Parkes 64-m telescope; and (3) the Effelsberg 100-m radio telescope.

Appendix A.1: GBT observations and data preparation

The GBT Hi maps are created from Hi spectral observations with the Green Bank Telescope, the details of which can be found in Boothroyd et al. (2011) and Martin et al. (in prep.). Over 800 deg² were mapped between 2002 and 2010, each following roughly the same observing strategy. The high Galactic latitude fields were mapped using scans along lines of constant declination or constant Galactic latitude at a scan rate equalling 3.5′ every 4 s. Each subsequent scan is offset by 3.5′ in the corresponding orthogonal direction. This strategy results in a rectangular region of the Hi sky sampled every 3.5′. Some of the regions have been scanned in this way multiple times, in order to increase the S/N, and to investigate the stability of the system (Boothroyd et al. 2011).

The spectra are first converted from their antenna temperature scale to a brightness temperature (T_b) scale. This involves calibration and stray radiation corrections, as discussed in Boothroyd et al. (2011). All T_b spectra for a corresponding region are assigned to a 3.5′ Sanson-Flamsteed-projection grid (SFL-projection) using convolution with an optimized tapered Bessel function, in order to minimize noise on spatial scales smaller than the beam (Mangum et al. 2007). Note that this observing strategy and gridding choice results in a final cube resolution of 9.55′×9.24′, slightly broader than the inherent 9.1′×9.0′ GBT 21-cm beam.

Each spectrum is recorded using in-band frequency-switching, resulting in velocity coverage of −450 km s^-1 ≤ V_LSR ≤ 355 km s^-1, with very flat baselines. Any residual baseline is removed on a pixel-by-pixel basis by fitting the emission-free channels in the final calibrated cube using a third-order polynomial.

The archival LH2 field was observed using a 3′ grid pattern with the GBT spectral processor and was calibrated accordingly. As the residual baseline behaviour is different for the archival data, only a linear polynomial was fit to the emission-free channels. The rest of the processing was identical to that described above.

The individual 0.8 km s^-1 spectral channels of the cubes are integrated to convert the T_b spectra into N_HI maps: ${\mathit{N}}_{\mathrm{HI}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\mathrm{1.823}\mathrm{\times }{\mathrm{10}}^{\mathrm{18}}\sum _{\mathit{v}}{\mathit{T}}_{\mathrm{b}}\mathrm{\left(}\mathit{x,y,v}\mathrm{\right)}\mathit{\tau }\mathrm{\left(}{\mathit{T}}_{\mathrm{s}}\mathrm{\right)}\mathit{\delta v,}$(A.1)where the sum is over a given velocity range, v, and δv is the 0.80 km s^-1 channel spacing. The quantity τ is the opacity correction for spin temperature, T_s: $\mathit{\tau }\mathrm{\left(}{\mathit{T}}_{\mathrm{s}}\mathrm{\right)}\mathrm{=}\mathrm{-}\ln \mathrm{(}\mathrm{1}\mathrm{-}{\mathit{T}}_{\mathrm{b}}\mathit{/}{\mathit{T}}_{\mathrm{s}}\mathrm{)}\mathit{.}$(A.2)For the adopted value of T_s = 80 K, these corrections are all less than 5% for our CIB fields (Planck Collaboration 2011RPlanck Collaboration XVIII 2011).

The velocity ranges over which the integrations are performed are selected using the observed velocity structure in each of the cubes. The models presented here subdivide each cube into three velocity-selected components: a local component; intermediate-velocity clouds, IVCs; and high-velocity clouds, HVCs. Divisions between components are distinguishable by reductions in structure (as measured through the standard deviation of individual channel maps) as one progresses through the data cube, channel by channel. More details can be found in Planck Collaboration 2011RPlanck Collaboration XVIII (2011).

Appendix A.2: GASS observations and data preparation

The GASS survey is a 21-cm line survey covering the southern sky for all declinations δ ≲ 1°. The observations were made with the multibeam system on the 64-m Parkes radio telescope. The intrinsic angular resolution of the data is 14.4′ full width at half maximum (FWHM). The velocity resolution is 1.0 km s^-1 and the useful bandpass covers a velocity range | v_lsr | ≲ 468 km s^-1 for all of the observations; some data cover up to | v_lsr | ≲ 500 km s^-1. GASS is the most sensitive, highest angular resolution large-scale survey of Galactic Hi emission ever made in the southern sky. The observations are described in McClure-Griffiths et al. (2009). We used data from the final data release (Kalberla et al. 2010) that were corrected for instrumental effects and radio-frequency interference (RFI). The data were gridded on a Cartesian grid on the Magellanic stream (MS) coordinate system as defined by Nidever et al. (2008). To minimize the noise and eliminate residual instrumental problems, we calculated a second 3D data cube with a beam of 0.5° FWHM, smoothing at the same time in velocity by 8 km s^-1. Emission below a 5σ level of 30 mK in the smoothed data-cube was considered as insignificant and was accordingly zeroed. For details in data processing and analysis see Venzmer et al. (2012).

When looking at the Hi data cube in the southern sky, one of the most prominent structures is linked to hydrogen gas in the Magellanic stream, in the disks of the Magellanic clouds, and in the stream’s leading arm. In particular, the Magellanic stream, stretches over 100° behind the Large and Small Magellanic clouds. We thus need to remove this contamination to be able to use the Hi as a tracer of Galactic dust.

Aiming to separate Galactic emission from the observations of the Magellanic system we calculated the expected Milky Way emission according to the model of Kalberla & Dedes (2008). Velocities for components in direction towards the southern Galactic pole were shifted by −5 km s^-1 to mimic the apparent infall (see Weaver 1974). Comparisons between the emission and the model at two particular velocities in the Galactic standard of rest frame are shown in Fig. A.1. The strong Galactic emission can be traced to weak extended line wings that are well represented by the model. It is therefore feasible to use the model to predict regions in the 3D data cube that are most probably occupied by Milky Way emission. We used a clip level of 60 mK, the lowest isophote in Fig. A.1 that delineates the disk emission. Such a treatment extracts most of the Galactic emission, however we must take into account the fact that at positions with strongly blended lines Magellanic emission also gets included. A higher clip level would minimize this problem, although, at the same time the wings of the Galactic emission would be affected. The chosen clip level of 60 mK (comparable to the instrumental noise), is a good compromise.

As the final step in the reduction of the GASS data we integrated the Hi emission over the appropriate velocity range for each individual position in the N_side = 512HEALPix database to obtain column densities of the Galactic gas. To avoid any interpolation errors, we extracted profiles from the original GASS database; all intermediate data products that have been described above served only to discriminate Galactic emission from the Magellanic stream. Due to the finite beam size of the Parkes telescope and the Gaussian weighting that was used for the gridding, the effective resolution of the HEALPix data is 16.2′ FWHM. Present a detailed analysis of the dust-Hi correlation in this field.

Appendix A.3: EBHIS observations and data analyses

The Effelsberg-Bonn Hi Survey (EBHIS) comprises an Hi survey of the entire northern sky for all declinations δ ≳ − 5° with the Effelsberg 100-m telescope. The bandwidth of 100 MHz covers − 1000 km s^-1 ≤ v_LSR ≤ 19 000 km s^-1. This allows us to study the detailed Milky Way Hi structure as well as the local Universe up to a redshift of z ≃ 0.07 (Kerp et al. 2011), with an effective velocity resolution of about 2.1 km s^-1.

We selected from the early survey data a clean high Galactic latitude field. The observations of the field of interest were performed during the summer of 2011. Following the standard observing strategy (Kerp et al. 2011) individual fields of 25 deg² were measured. In addition to the data reduction and calibration pipeline of Winkel et al. (2010), the Hi data were corrected with an improved RFI mitigation detection algorithm and the absolute offsets between the individual fields were minimized to a level of N_HI ≤ 3 × 10¹⁸ cm^-2. As for the GBT fields, the EBHIS data were put on a 3.5′ SFL-projection grid.

The overall Hi emission in this field can be characterized by intermediate- and low-velocity emission populating the radial velocity range − 80 km s^-1 ≤ v_LSR ≤ + 20 km s^-1.

The area of interest is about 130 deg² with Hi column densities below N_HI = 2 × 10²⁰ cm^-2. The total Hi column density range is 0.98 × 10²⁰ cm^-2 ≤ N_HI ≤ 7.5 × 10²⁰ cm^-2. Accordingly we expect infrared excess emission associated with molecular hydrogen towards areas of N_HI ≥ 3 × 10²⁰ cm^-2, while the low-column density regions should allow us to study the emission associated with the CIB.

In the first step we performed a linear fit of the EBHIS Hi data integrated across the velocity range − 160 km s^-1 ≤ v_LSR ≤ + 160 km s^-1 to the HFI 857 GHz map: I₈₅₇ = a × N_HI + b, for N_HI ≤ 2 × 10²⁰ cm^-2. The residual map shows dust excess, as expected when the total infrared emissivity traces Hi and H₂ (Planck Collaboration 2011SPlanck Collaboration XIX 2011; Planck Collaboration 2011XPlanck Collaboration XXIV 2011). In agreement with our expectation dust excess clouds show up with Hi column densities of N_HI ≥ 3 × 10²⁰ cm^-2. But surprisingly, the residual map shows dust-deficient Hi clouds that populate the column density range 2.0 ≤ N_HI ≤ 5.5 × 10²⁰ cm^-2. The Hi emission of these dust-deficient clouds can be characterized by brightness temperatures T_B> 15 K and Δν(FWHM) ≤ 5 km s^-1. These Hi lines are associated with individual, cold IVCs. A positional correlation between the dust-deficient regions and the IVCs show a very strong correlation. We evaluated the dust deficiency of the IVCs and found that they are not “under-luminous”, but infrared dark.

To extract the CIB in this field we produced a mask excluding the excess as well as the dust-deficient clouds. About 90 deg² are suitable for CIB analysis.

Appendix B: From infrared luminosity density to emissivities

The emissivity, often written j_ν, is given by (Pénin et al. 2012a) ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{d}\mathit{\chi }\mathit{/}\mathrm{d}\mathit{z}}\mathrm{\int }{\mathit{S}}_{\mathit{\nu }}\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{N}}{\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{z}}\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\mathit{,}$(B.1)where χ is the comoving distance, S_ν the flux density, and d²N/dS_ν dz the number of sources per flux density and redshift interval. If there are several types of galaxy (labelled by “t”) with different SEDs⁹, one can rewrite Eq. (B.1) as ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{d}\mathit{\chi }\mathit{/}\mathrm{d}\mathit{z}}\sum _{\mathit{t}}\mathrm{\int }{\mathit{S}}_{\mathit{\nu }}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{t}}}{\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{z}}\mathrm{d}{\mathit{S}}_{\mathit{\nu }}\mathit{.}$(B.2)We then introduce s_ν,t the flux density of an L_IR = 1 L_⊙ source with an SED of a given type. This quantity varies only with z for a given type of SED. We can thus easily change the variable S_ν to L_IR in the integral: ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{d}\mathit{\chi }\mathit{/}\mathrm{d}\mathit{z}}\sum _{\mathit{t}}{\mathit{s}}_{\mathit{\nu }}\mathrm{\int }{\mathit{L}}_{\mathrm{IR}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{t}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{z}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{.}$(B.3)We can slightly modify this expression so that it contains the bolometric infrared luminosity function of each galaxy type d²N_t/ dL_IR dV: ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}\mathrm{+}\mathit{z}}{\mathrm{d}\mathit{\chi }\mathit{/}\mathrm{d}\mathit{z}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}}\sum _{\mathrm{t}}{\mathit{s}}_{\mathit{\nu }}\mathrm{\int }{\mathit{L}}_{\mathrm{IR}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{t}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{V}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{,}$(B.4)where V is the comoving volume. We can simplify this expression assuming a flat ΛCDM Universe with $\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathit{c}}{{\mathit{H}}_{\mathrm{0}}\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}$(B.5)and $\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathit{c}{\mathit{\chi }}^{\mathrm{2}}}{{\mathit{H}}_{\mathrm{0}}\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\Lambda }}}}\mathrm{=}\frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}\mathit{z}}{\mathit{\chi }}^{\mathrm{2}}\mathit{.}$(B.6)We can thus simplify Eq. (B.4): ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{\chi }}^{\mathrm{2}}\sum _{\mathit{t}}{\mathit{s}}_{\mathit{\nu }}\mathrm{\int }{\mathit{L}}_{\mathrm{IR}}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{t}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{V}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{.}$(B.7)We introduce the contribution to the infrared luminosity density of a given type of SED: ${\mathit{\rho }}_{\mathrm{IR}\mathit{,t}}\mathrm{=}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{N}}_{\mathit{t}}}{\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\hspace{0.17em}\mathrm{d}\mathit{V}}\mathrm{d}{\mathit{L}}_{\mathrm{IR}}\mathit{,}$(B.8)where ρ_IR is the total infrared luminosity density. We can assume a simple conversion between L_IR and ρ_SFR, the star formation rate density, using the Kennicutt (1998) constant K. Then Eq. (B.7) can be rewritten as: ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{\chi }}^{\mathrm{2}}\frac{{\mathit{\rho }}_{\mathrm{SFR}}}{\mathit{K}}\mathrm{\times }\left(\frac{{\sum }_{\mathit{t}}{\mathit{s}}_{\mathit{\nu }}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{IR}\mathit{,t}}}{{\sum }_{\mathit{t}}{\mathit{\rho }}_{\mathrm{IR}\mathit{,t}}}\right)\mathit{,}$(B.9)where ∑ _ts_νρ_IR,t/ ∑ _tρ_IR,t is the effective SED of infrared galaxies (noted s_ν,eff), i.e., the mean SED using a weighting by the contribution to the infrared luminosity density. Here we use the model of Béthermin et al. (2012a) to compute the contribution at each redshift. This model is based on the observed strong correlation between stellar mass and star formation rate in these galaxies (called main-sequence galaxies), and new SED templates from Magdis et al. (2012). We finally obtain: ${\mathit{j}}_{\mathit{\nu }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{SFR}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\hspace{0.17em}{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\hspace{0.17em}{\mathit{s}}_{\mathit{\nu ,}\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{K}}\mathrm{·}$(B.10)

Appendix C: Effective bias notion for the linear model

We discuss here the link between the approach chosen to model the 2-halo term in Sects. 5.4 (linear model) and 5.5 (HOD model). We consider for simplicity only the case of the auto-spectrum, but this can be easily generalized to the cross-spectrum. Equation (48) of the HOD model can be re-written as ${\mathit{P}}_{\mathrm{2}\mathrm{h}}^{\mathit{\nu \nu }}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{=}\mathrm{\int }{\left(\frac{\mathit{a}}{\mathit{\chi }}\right)}^{\mathrm{2}}{\left(\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\right)}^{\mathrm{2}}{\mathit{P}}_{\mathrm{lin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{d}\mathit{\chi ,}$(C.1)where dj/dM is the differential contribution of haloes of mass M, at redshift z, to the emissivity. It can be computed using $\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{M}}\mathrm{\times }\mathrm{\left\{}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{cen}}\mathrm{\right(}\mathit{M,z}\mathrm{)}\mathrm{+}{\mathit{f}}_{\mathit{\nu }}^{\mathrm{sat}}\mathrm{(}\mathit{M,z}{\mathrm{\left)}}^{\mathrm{\right\}}}\mathit{.}$(C.2)The total emissivity is thus $\mathit{j}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{d}\mathit{M}\mathit{.}$(C.3)In order to simplify Eq. (C.1), we introduce an effective bias, ${\mathit{b}}_{\mathrm{eff}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\int }\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathit{b}\mathrm{\left(}\mathit{M,z}\mathrm{\right)}\mathrm{d}\mathit{M}}{\mathrm{\int }\frac{\mathrm{d}\mathit{j}}{\mathrm{d}\mathit{M}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}\mathrm{d}\mathit{M}}\mathrm{·}$(C.4)This effective bias is the mean bias of haloes hosting the infrared galaxies, weighted by their differential contribution to the emissivities. We finally obtain a simpler equation, which is used in the linear model: ${\mathit{P}}_{\mathrm{2}\mathrm{h}}^{\mathit{\nu \nu }}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{=}\mathrm{\int }{\left(\frac{\mathit{a}}{\mathit{\chi }}\right)}^{\mathrm{2}}{\mathit{j}}^{\mathrm{2}}\mathrm{\left(}\mathit{\nu ,z}\mathrm{\right)}{\mathit{b}}_{\mathrm{eff}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{lin}}\mathrm{\left(}\mathit{k,z}\mathrm{\right)}\mathrm{d}\mathit{\chi }\mathit{.}$(C.5)

Appendix D: Tables

All Tables

All Figures