© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration I 2014–Planck Collaboration XXXI 2014), describes the CMB power spectra and the corresponding likelihood that we derive from the Planck data.

The power spectrum of the cosmic microwave background (CMB) is a unique signature of the underlying cosmological model (e.g., Spergel et al. 2003; Hinshaw et al. 2013). It has been measured over the whole sky by COBE and WMAP, and over smaller regions by ground-based and sub-orbital experiments (e.g., Tristram et al. 2005; Jones et al. 2006; Reichardt et al. 2009; Fowler et al. 2010; Das et al. 2011; Keisler et al. 2011; Story et al. 2013; Das et al. 2014). By mapping the whole sky to scales of a few arc minutes, Planck now measures the power spectrum over an unprecedented range of scales from a single experiment. To estimate cosmological parameters from the power spectrum requires a likelihood function that propagates uncertainties.

In this paper we describe the power spectra obtained from the Planck temperature data, as well as the associated likelihood function. Since the probability distribution of the power spectrum is non-Gaussian at large scales, we follow a hybrid approach to construct this likelihood (Efstathiou 2004, 2006), using a Gibbs sampling based approach at low multipoles, ℓ, and a pseudo-C_ℓ technique at high multipoles (Hivon et al. 2002) as well as an approximation to the quadratic maximum likelihood estimator (Rocha et al. 2011, 2010) at high multipoles.

The high-ℓ part of the Planck likelihood (ℓ ≥ 50) is based on power spectra estimated from each Planck detector in the frequency range 100 to 217 GHz, allowing careful assessment of each detector’s response to the sky emission. We implement three independent likelihood methods. The first, used in the distributed likelihood code, estimates the power spectrum at every multipole, together with the associated covariance matrix. The second takes a simplified form, binning the spectra, and is used to explore the stability of the results with respect to different instrumental and astrophysical systematic effects, the third one takes a CMB map estimated from component separation techniques (hence using all Planck frequency maps) and estimates the band powers and associated covariance matrices. The methods give consistent results.

Unresolved extragalactic foregrounds make a significant contribution to the power spectra at high multipoles. We develop a model for these foregrounds, designed to allow the Planck likelihood to be combined with high resolution data from the Atacama Cosmology Telescope (ACT) and the South Pole Telescope (SPT). We combine frequencies and model unresolved foregrounds in a physical way, as in e.g., Shirokoff et al. (2011); Dunkley et al. (2011); Reichardt et al. (2012), performing component separation at small scales at the power spectrum level. On large scales, ℓ < 50, Galactic contamination is more significant. We use the Planck temperature maps in the range 30 ≤ ν ≤ 353 GHz to separate Galactic foregrounds in the maps, and then estimate the full probability distribution of the CMB power spectrum.

This paper is structured as follows. In Sect. 2 we describe the pseudo-C_ℓ likelihoods, and in Sect. 3 set up the foreground model. The power spectra and derived cosmological parameters are presented in Sects. 4 and 5, and an assessment of their accuracy and robustness is made in Sects. 6 and 7. In Sect. 8 we describe the low-ℓ likelihood, and conclude by presenting the complete Planck likelihood in Sect. 9.

2. High-ℓ likelihoods

The Planck maps consist of the order 5 × 10⁷ pixels for each detector (Górski et al. 2005), so a likelihood described directly at the pixel level would be too time consuming. A significant compression of data can be achieved with minimal information loss using pseudo-C_ℓ power spectra, even in the case of incomplete sky coverage. Here we describe the form of the likelihood function of the compressed data, given a sky signal and instrumental model.

Following Hamimeche & Lewis (2008), we assume a Gaussian form of the likelihood (Bond et al. 2000) based on pseudo-spectra that have been corrected to account for partial sky masking (Hivon et al. 2002) (for an overview on high-ℓ Likelihoods see e.g. Rocha et al. 2010). We use a “fine-grained” data description, computing spectra of maps from individual detectors or detector sets. Table 1 describes the 13 maps used in the analysis, spanning 100 to 217 GHz, which are compared with the general purpose “ds0” maps used in input to the component separation methods. These maps are based on the Planck nominal mission. In the case of the polarised detector sets (PSB type in Table 1), the map-making algorithm uses the four combined detectors to produce maps of the T, Q, U Stokes parameters, but only the temperature maps produced from each PSB detector set are used in the likelihood analyses. We compute the spectra at these multiple frequencies to simultaneously constrain the CMB and foreground contributions. We choose these frequencies as a trade-off between adding further information, and adding further complexity to the foreground model, which would be needed to include the adjacent 70 and 353 GHz channels (see Sect. 3 for further discussion). In our baseline analysis the spectra are computed at each multipole, together with an estimate of the full covariance matrix with off-diagonal errors between different spectra and multipoles. As in the WMAP analysis, we use only cross-spectra between detectors, alleviating the need to accurately model the mean noise contribution.

In this section, we begin with a reminder of the pseudo-spectrum approach, and describe our baseline likelihood distribution, hereafter referred to as the CamSpec likelihood. We then show how a compression of spectra within a given frequency can be achieved with negligible loss of information. We describe the signal and instrument model, including detector noise properties, calibration, and beam uncertainties.

Next, we describe an alternative, simpler, form of the likelihood, hereafter referred to as Plik, based on binned power spectra with an inverse-Wishart distribution. This does not require the pre-computation of large covariance matrices, so changing the sky or instrument modelling is straightforward. This simpler form of the likelihood is used to assess the robustness of our likelihood methodology with respect to technical choices and astrophysical foreground modelling.

In Sect. 7.5 we also compare these likelihoods to the Likelihood estimated from band power spectra computed directly from CMB maps estimated by multi-frequency component separation (Planck Collaboration XII 2014). The band powers and covariance matrices are estimated with the XFaster approach (Rocha et al. 2011, 2010).

2.1. The CamSpec likelihood

The temperature maps produced from individual detector sets are first weighted with apodised Galactic and point source masks that are common within a given frequency band (see Sect. 3) for details. On top of these masks crafted to reject astrophysical contamination, each detector set has a small number of unobserved pixels that are added to the corresponding detector set masks.

Let us define $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathit{\ell m}}^{\mathit{i}}\end{array}$ as the spherical harmonic coefficients of the weighted temperature map of detector set i. The pseudo-spectrum at multipole ℓ, for the detector pair (i,j), is then given by $\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{\mathit{ij}}}_{\mathrm{˜}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}}\stackrel{}{{{\mathit{T}}_{\mathit{\ell m}}^{\mathit{i}}}_{\mathrm{˜}}}\stackrel{}{{{\mathit{T}}_{\mathit{\ell m}}^{\mathit{j}†}}_{\mathrm{˜}}}\mathit{,}$(1)where the dagger, †, denotes the Hermitian transpose. This is related to the “deconvolved” spectrum, Ĉ^T_ij, by a coupling matrix, $\stackrel{}{{{\mathit{C}}^{{\mathit{T}}_{\mathit{ij}}}}_{\mathrm{˜}}}\mathrm{=}{\mathit{M}}_{\mathit{ij}}^{\mathit{TT}}\mathit{Ĉ}{\mathit{T}}_{\mathit{ij}}\mathit{.}$(2)For an isotropic signal on the sky, the ensemble average of these deconvolved spectra are equal to the spectra of the theoretical models (including CMB and isotropic unresolved foregrounds) that we wish to test. For completeness, the coupling matrices are given explicitly in Appendix A.1.

In the first method, CamSpec, we form the deconvolved spectra Ĉ_ℓ without any prior smoothing of the pseudo-spectra $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}\end{array}$. Even for the largest sky masks used in our analysis (see Sect. 3), the coupling matrices are non-singular. The deconvolution requires the evaluation of ~${\mathit{N}}_{\mathrm{map}}^{\mathrm{2}}$ coupling matrices for a data set with N_map sky maps, which takes a moderate, but not excessive, amount of computer time.

A more challenging computational task is to compute the covariances of the pseudo-spectra, i.e., $\mathrm{Cov}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}^{{\mathit{T}}_{\mathit{ij}}}\end{array}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{C}}^{{\mathit{T}}_{\mathit{pq}}}\end{array}\mathrm{\right)}\mathit{.}$ Here we need to compute ${\mathit{N}}_{\mathrm{map}}^{\mathrm{4}}$ coupling matrices, and the problem rapidly becomes computationally intractable even for relatively low values of N_map. For the moment we assume that these covariance matrices are available and describe their computation in Appendix A.4. We use the notation $\begin{array}{c}{}_{{˜}}\\ {X}\end{array}\mathrm{=}\mathrm{Vec}\mathrm{\left(}\begin{array}{c}{}_{{˜}}\\ {C}\end{array}\mathrm{\right)}$ to denote a column vector for which the index p of a single element X_p denotes the map combination (i,j) and multipole ℓ. We denote the covariance matrix of this vector as $\begin{array}{c}\mathrm{˜}\\ \mathrm{ℳ}\end{array}\mathrm{=}\mathrm{(}\begin{array}{c}{˜}\\ {X}\end{array}\mathrm{-}\mathrm{⟨}\begin{array}{c}{˜}\\ {X}\end{array}\mathrm{⟩}\mathrm{)}\mathrm{(}\begin{array}{c}{˜}\\ {X}\end{array}\mathrm{-}\mathrm{⟨}\begin{array}{c}{˜}\\ {X}\end{array}\mathrm{⟩}{\mathrm{)}}^{\mathrm{T}}\mathit{.}$(3)As explained later, the deconvolved detector set cross-spectra given by Eq. (2) can be efficiently combined within a given frequency pair after a small effective recalibration, taking their respective isotropised effective beam transfer function and noise levels into account (see Appendix A.3 for the detailed procedure). Here, like in all Planck papers, the effective beam means the beam, which accounts for the combined effect of the instrument and of the complete processing chain. Covariance estimates of these combined spectra can be deduced from those of the detector set cross-spectra. The covariance matrix is computed for a fixed fiducial model, and we approximate the likelihood as a Gaussian, described in Appendix A.5. The likelihood thus takes the form p = e^{− S} with $\begin{array}{ccc}\mathit{S}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{(}{X̂}\mathrm{-}{{X}}^{\mathrm{)}}}^{\mathrm{T}}\widehat{\mathrm{ℳ}}-1\mathrm{(}{X̂}\mathrm{-}{{X}}^{\mathrm{)}}\mathit{.}& & \end{array}$For the current analysis we include the following (deconvolved) spectrum combinations, ${X̂}\mathrm{=}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}\\ \mathit{\ell }\end{array}\mathit{,}\mathit{Ĉ}\begin{array}{c}\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{143}\\ \mathit{\ell }\end{array}\mathit{,}\mathit{Ĉ}\begin{array}{c}\mathrm{217}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}\\ \mathit{\ell }\end{array}\mathit{,}\mathit{Ĉ}\begin{array}{c}\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}\\ \mathit{\ell }\end{array}\mathrm{\right)}\mathit{,}$(4)coupled to a parametric model of the CMB and foreground power spectra. The multipole ranges we select, as well as the masks, depend on frequency, as described in Sect. 5. The tailored selection of multipole ranges and masks we adopt here is designed to keep the Galactic contamination as low as possible in both real and harmonic space, while keeping a reasonably low cosmic variance. This procedure allows us to keep the Galactic emission model simple (see Sect. 3.2). We do not include the 100 × 143 and 100 × 217 spectra since these spectra carry little additional information about the primary CMB anisotropies, but would require us to solve for additional unresolved foreground parameters. This tradeoff of information versus complexity was also considered for the use of the 70 GHz and 353 GHz data, which we choose not to include except for cross-checks.

The fiducial covariance matrix is composed of the blocks shown in Fig. 1. The off-diagonal blocks in this matrix accurately account for the correlations between the power spectra at different frequencies.

This description would be sufficient for perfectly known calibrations and beam transfer functions of each detector sets’ cross-spectra. Planck Collaboration VII (2014) describes in detail these uncertainties, and shows that for each detector set pair, (i,j), the effective beam transfer function can be expressed as ${\mathit{B}}^{\mathit{ij}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{mean}}^{\mathit{ij}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\exp \left(\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{modes}}}{\mathit{g}}_{\mathit{k}}^{\mathit{ij}}{\mathit{E}}_{\mathit{k}}^{\mathit{ij}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right)\mathit{,}$(5)described further in Appendix A.6, with n_modes beam error eigenmodes ${\mathit{E}}_{\mathit{k}}^{\mathit{ij}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$, and their covariance matrix. These modes are then combined into generalised beam eigenmodes corresponding to the spectra . The associated covariance matrix is used to construct a Gaussian posterior distribution of the eigenmodes, which allows marginalisation over the beam uncertainties.

Finally, in the construction of the covariance matrix, one needs to accurately specify the contribution of the instrumental noise. Even if there is no bias on the spectra due to instrumental noise, having removed auto-spectra, the latter dominates the covariance matrix on small scales. Fortunately, the Planck scanning strategy at the ring level allows us to make estimates of the noise pseudo-spectra from half-ring difference maps (see Planck Collaboration XII 2014). These half-ring difference maps, together with the knowledge of the noise variance per pixel for each detector set, can be used to derive the noise contribution to the covariance matrix with good accuracy (see Appendix A.8 for details).

2.2. The Plik likelihood

We now describe the alternative form of the likelihood, inspired by Cardoso et al. (2008), used for cross-checks and robustness tests. We start from the full-sky exact likelihood for a Gaussian signal, which for N_map detector maps is given by $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathrm{maps}\mathrm{\right|}\mathit{\theta }\mathrm{)}\mathrm{\propto }\exp \mathrm{-}\left\{\sum _{\mathit{\ell }}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{𝒦}(\mathit{Ĉ}\mathit{\ell }\mathit{,}\hspace{0.17em}{{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)})\right\}\mathit{,}& & \end{array}$where θ is a vector containing the parameters of the signal model, and Ĉ_ℓ are the empirical angular spectra. denotes the Kullback divergence between two n-variate zero-mean Gaussian distributions with covariance matrices A and B, and is given by $\begin{array}{ccc}\mathrm{𝒦}({A}\mathit{,}\hspace{0.17em}{B})\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}[\mathrm{tr}({A}{{B}}^{-1})\mathrm{-}\log \begin{array}{c}\det \end{array}({A}{{B}}^{-1})\mathrm{-}\mathit{n}]\mathit{.}& & \end{array}$As already noted, a sky cut introduces off-diagonal couplings between different multipoles. In this method we bin the power spectra in such a way that these off-diagonal terms of the covariance are negligible. This is adequate to model sources with slowly varying spectra, such as foregrounds, and the CMB anisotropies for standard cosmologies. In this case, the likelihood takes the form $\begin{array}{ccc}\mathit{p}\mathrm{\left(}\mathrm{maps}\mathrm{\right|}\mathit{\theta }\mathrm{)}\mathrm{\propto }\exp \mathrm{-}\mathrm{ℒ}\mathrm{(}\mathit{\theta }\mathrm{\left)}\mathit{,} \mathrm{with} \mathrm{ℒ}\mathrm{\right(}\mathit{\theta }\mathrm{)}\mathrm{=}\sum _{\mathit{q}\mathrm{=}\mathrm{1}}^{\mathit{Q}}{\mathit{n}}_{\mathit{q}}\hspace{0.17em}\mathrm{𝒦}(\mathit{Ĉ}\mathit{q}\mathit{,}\hspace{0.17em}{{C}}_{\mathit{q}})\mathit{,}& & \end{array}$(6)where the angular spectra Ĉ_ℓ for each cross-frequency spectrum have been averaged into Q spectral bins using spectral windows w_q(ℓ) (q = 1,...,Q), with $\begin{array}{ccc}\widehat{{\mathit{C}}_{\mathit{q}}}\mathrm{=}\sum _{\mathit{\ell }}{\mathit{w}}_{\mathit{q}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{}_{\mathit{\ell }}\widehat{{\mathit{C}}_{\mathit{\ell }}}\mathit{,} {{C}}_{\mathit{q}}\mathrm{=}\sum _{\mathit{\ell }}{\mathit{w}}_{\mathit{q}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{{C}}_{\mathit{\ell }}\mathit{.}& & \end{array}$Here w_q(ℓ) denotes the window function for the qth bin, and the same symbol, C, is used to denote binned or unbinned spectra. The effective number of modes in the qth bin is $\begin{array}{ccc}{\mathit{n}}_{\mathit{q}}\mathrm{=}{\mathit{f}}_{\mathrm{sky}}\mathrm{·}\frac{\mathrm{\left(}{\sum }_{\mathit{\ell }}{\mathit{w}}_{\mathit{q}}\mathrm{\right(}\mathit{\ell }{\mathrm{)}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{{\sum }_{\mathit{\ell }}{\mathit{w}}_{\mathit{q}}\mathrm{(}\mathit{\ell }{\mathrm{)}}^{\mathrm{4}}\mathit{/}\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{·}& & \end{array}$We adopt a spectral binning defined by $\begin{array}{ccc}{\mathit{w}}_{\mathit{q}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}& & \end{array}$The Plik bin width is Δℓ = 9 from ℓ = 100 to ℓ = 1503, then Δℓ = 17 to ℓ = 2013, and finally Δℓ = 33 to ℓ_max = 2508. This ensures that correlations between any two bins are smaller than 10%.

While this binned likelihood approximation does not fully capture all couplings between different multipoles, it has a notable advantage in computational speed, and it agrees well with the primary likelihood. It is therefore very well suited for performing an extensive suite of robustness tests, as many more parameters can be considered in a short time. Further, instrumental effects can be investigated quickly to assess the agreement between pairs of detectors within a frequency channel, such as individual detector calibrations and beam errors.

A specific example is the impact of (effective) beam uncertainty parameters on the likelihood. This can be investigated by re-expressing the model covariance matrices as ${{C}}_{\mathit{\ell }}\mathrm{=}{{B}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\hspace{0.17em}{{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}{{B}}_{\mathit{\ell }}\mathrm{(}\mathit{\gamma }{\mathrm{)}}^{\mathrm{T}}\mathit{,}$(7)where C_ℓ(θ) is the model covariance including both signal and noise, and B_ℓ(γ) is a diagonal matrix encoding the beam and calibration errors with elements given by²${\mathit{B}}_{\mathit{\ell }}^{\mathit{i}}\mathrm{\left(}\mathit{\gamma }\mathrm{\right)}\mathrm{=}\exp \left(\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{modes}}}{\mathit{\delta }}_{\mathit{k}}^{\mathit{i}}\hspace{0.17em}{\mathit{E}}_{\mathit{k}}^{\mathit{ii}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\right)\mathit{.}$(8)Here, ${\mathit{E}}_{\mathit{k}}^{\mathit{ii}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$ are the eigenmodes of the (auto-)spectra, similar to Eq. (5). Note that Eq. (7) does not contain the mean beam transfer function, since it is already included in the empirical spectra. Thus, using Eq. (7) Plik approximates the cross-spectrum beam errors as the harmonic mean of the corresponding auto-spectrum beam errors, under the assumption that B_ℓ is diagonal between detectors. This approximate factorisation is intrinsically linked to the assumed Kullback shape of the Plik likelihood, and is later demonstrated to work well for both simulations and data.

The Plik likelihood method also provides a direct estimate of the detector noise power spectra as it can include the empirical auto-spectra, and we find that these noise estimates are in good agreement with the noise spectra used to construct the CamSpec likelihood covariance matrix. The method can also produce a binned CMB power spectrum independent of the underlying cosmological model, providing a direct quality assessment of the foreground model parametrisation. In practice, we proceed in two steps. First, we jointly estimate the noise together with all other parameters using both auto and cross-spectra. Then we fix the noise estimates, and use the fiducial Gaussian approximation to explore the remaining free parameters excluding the auto-spectra, optionally including only specific data combinations.

3. Foreground emission model and sky masks

3.1. Sky masks

The Galactic emission varies strongly in both complexity and strength across the sky. It is therefore necessary to find a balance between maximizing the sky coverage to reduce statistical uncertainties, and establishing a simple yet efficient foreground model. In this paper, we threshold an ILC-subtracted 353  GHz temperature map to define a basic set of diffuse Galactic masks, which form a sequence of increasing sky fraction, to minimise the contribution from diffuse dust emission. We refer to them using the percentage of the sky retained: G22, G35, G45, G56, G65. We also use a point source mask, labelled PS96, which is based on the union of the point sources detected from the channels in the range 100 to 353 GHz. The point source flux cut is not critical, since the amplitudes of the Poisson contributions of unresolved sources are allowed to vary over a wide range in the likelihood analysis. Thus, we do not impose tight priors from source counts and other CMB experiments on the Poisson amplitudes. All theses masks are shown in Fig. B.1.

For Planck, we need to estimate the covariance matrices to percent level precision. For temperature spectra, and in the absence of point source holes, this precision can be achieved with sharp, non-apodised Galactic masks (Efstathiou 2004). However, the inclusion of point source holes introduces non-negligible low-ℓ power leakage, which in turn can generate errors of several percent in the covariance matrices. In order to avoid such power leakage, we also derive a series of apodised masks. For the Galactic masks, we proceed as follows. First, we smooth each mask with a five-degree Gaussian beam, and zero any pixels below a threshold of 0.15 (chosen to approximately keep the zeroed area the same as that of the original, unapodised mask). We then subtract 0.15 from the remaining pixels, and rescale the resulting mask by 1 / (1 − 0.15). This procedure is iterated a further two times. In order to retain sufficient sky area for the most conservative sky mask, a slightly less aggressive version of mask G22 was used to seed the apodisation process for that case. Each point source is apodised to 30′ FWHM, resulting in the PSA82 point source mask. The resulting set of apodised masks are shown in Fig. B.2.

For all cosmology analyses, we use three of the union of one of the apodised galactic masks and of the apodised point sources mask – CL31, CL39, and CL49, which are shown in Fig. 2. Table 2 summarizes the various masks and the corresponding sky fraction they retain.

3.2. Galactic emission

The contamination from diffuse Galactic emission at low to intermediate multipoles can be reduced to low levels compared to CMB anisotropies by a suitable choice of masking. However, even with conservative masking, the remaining Galactic emission at high multipoles is non-negligible compared to other unresolved components, such as the cosmic infrared background (CIB) anisotropies at 143 and 217 GHz. A clear way of demonstrating this is by differencing the power spectra computed with different masks, thereby highlighting the differences between the isotropic and non-isotropic unresolved components.

Let us define , a notation we shall use throughout this paper. Figure 3 shows (up to ℓ ≤ 1400) the 217 GHz power spectrum difference for the union masks mask0 ≡ G22∪PSA82 and mask1 ≡ G35∪PSA82, minus the corresponding difference for the 143 GHz frequency channel. Any isotropic contribution to the power spectrum (CMB, unresolved extragalactic sources etc.) cancels in such a mask-differenced power spectrum, leaving a non-isotropic signal of Galactic origin. Nevertheless, at low multipoles, a mask differenced power spectrum is still severely affected by cosmic variance of the primordial CMB; but this can be eliminated by forming a “double-difference” power spectrum between frequencies. This is what is plotted at ℓ ≤ 1400 in Fig. 3, rescaling to account for the small dust contribution at 143 GHz. Above ℓ = 1400, the noise in the 143 GHz power spectra becomes significant, and so we plot directly the 217 mask-differenced power spectrum. The figure gives an accurate estimate of the dust contribution within mask1 to the 217 GHz spectra over the entire multipole range plotted. In the same figure, these difference spectra are compared to the unbinned mask-differenced 857 GHz power spectrum, scaled to 217 GHz adopting a multiplicative factor³ of (9.93 × 10^-5)²; the dotted line shows a smooth fit to the unbinned spectrum. The agreement between this prediction and the actual dust emission at 217 GHz is excellent, and this demonstrates conclusively the existence of a small-scale dust emission component with an amplitude of ~5 − 15 μK² at 217 GHz if mask1 is used.

For cosmological parameter analysis this small-scale dust component must be taken into account, and several approaches may be considered:

• 1.

  Fit to a template shape, e.g., as shown by the dotted line inFig. 3.

• 2.

  Reduce the amplitude by further masking of the sky.

• 3.

  Attempt a component separation by using higher frequencies.

The main disadvantage of the third approach is a potential signal-to-noise penalty, depending on which frequencies are used, as well as confusion with other unresolved foregrounds. This is particularly problematic with regards to the CIB, which has a spectrum very similar to that of Galactic dust. In the following we therefore adopt the two former solutions.

It is important to understand the nature of the small scale dust emission, and, as far as possible, to disentangle this emission from the CIB contribution at the HFI cosmological frequencies. We use the 857 GHz power spectrum for this purpose, noting that the dust emission at 857 GHz is so intense that this particular map provides an effectively noise-free dust emission map. In Fig. 4 we again show the 857 GHz mask power spectrum difference, but this time plotted on a log-log scale. The solid line shows the corresponding best-fit model defined by ${\mathrm{𝒟}}_{\mathit{\ell }}\mathrm{=}\frac{\mathit{A}\hspace{0.17em}\mathrm{(}\mathrm{100}\mathit{/}\mathit{\ell }{\mathrm{)}}^{\mathit{\alpha }}}{\mathrm{[}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\ell }\mathit{/}{\mathit{\ell }}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{2}}{\mathrm{]}}^{\mathit{\gamma }\mathit{/}\mathrm{2}}}\mathit{,}$(9)with A = 5.729 × 10⁸μK², α = 0.169, ℓ_c = 90.5, and γ = 0.427. At high multipoles this fit asymptotically approaches C_ℓ ∝ ℓ^-2.6, which is compatible with previous knowledge about diffuse Galactic emission, i.e., a power-law behaviour with an index close to − 3 extending to high multipoles (see e.g., Miville-Deschênes et al. 2007). Note that throughout the paper, we use CMB thermodynamic units.

The upper panel in Fig. 5 shows the 857 GHz spectra for the four Galactic masks (G22, G35, G45 through to G56) with the point-source mask applied. They are compared to the 857 GHz CIB power spectrum from Planck Collaboration XVIII (2011), which for ℓ > 500 can be described approximately as ℓC_ℓ ≈ 7 ×10⁶μK². The best fit models to Eq. (9) are also shown (with α, ℓ_c and γ fixed to the values derived above, but with the amplitude allowed to vary) fitted to ℓ ≤ 500, where we expect diffuse dust emission to be dominant. The lower panel of Fig. 5 shows the same power spectra after subtracting the best-fit dust model. After subtracting the Galactic dust component, the recovered power spectra are consistent with the CIB measured in Planck Collaboration XVIII (2011) for all masks. The excess at high multipoles may be due to a combination of aliasing of large scale power through the point source masks at ℓ ≳ 3000, Galactic point sources, and uncertainties in the 857 GHz beams.

The model explains by construction the “double-difference” plot shown in Fig. 3. Specifically, this emission is consistent with a cirrus-like power spectrum, C_ℓ ∝ ℓ^-2.6, extrapolated to high multipoles. Furthermore, the results of Fig. 5 demonstrate that over a wide area of sky, we can understand the 857 GHz power spectrum in terms of a “universal” cirrus spectrum together with an isotropic CIB component. These results provide strong evidence that an extragalactic CIB component dominates over the diffuse Galactic emission at multipoles ℓ > 700 over the full range of HFI frequencies outside the G22 region (i.e. over the “cleanest” ~20% of sky). Although this is demonstrated for 857 GHz in Fig. 5, it applies over the entire HFI frequency range since Galactic dust and CIB have similar spectra. We therefore use the dust template of of Eq. (9) with suitably adjusted scaling factors to remove the high-ℓ contribution from dust from the empirical 143 × 143, 143 × 217 and 217 × 217 spectra, which are used in CamSpec.

We take a different approach for the Galactic dust correction with the Plik likelihood. Rather than correcting the empirical spectra during a pre-processing step, the Plik likelihood implements an explicit one-parameter model that describes the dust contribution to the cross-spectrum between detectors i and j, ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{Dust}}\mathrm{(}\mathit{i,}\hspace{0.17em}\mathit{j}\mathrm{)}\mathrm{=}{\mathit{A}}^{\mathrm{Dust}}\hspace{0.17em}\mathit{F}\mathrm{\left(}{\mathit{\nu }}_{\mathit{i}}\mathit{,}{\mathit{\nu }}_{\mathrm{0}}\mathrm{\right)}\mathit{F}\mathrm{\left(}{\mathit{\nu }}_{\mathit{j,}}{\mathit{\nu }}_{\mathrm{0}}\mathrm{\right)}\hspace{0.17em}{\left(\frac{\mathit{\ell }}{\mathrm{500}}\right)}^{\mathrm{-}{\mathit{\gamma }}_{\mathrm{d}}}\hspace{0.17em}{\mathit{g}}_{\mathit{i}}^{\mathrm{Dust}}{\mathit{g}}_{\mathit{j}}^{\mathrm{Dust}}\mathit{.}$(10)Here $\mathit{F}\mathrm{\left(}{\mathit{\nu }}_{\mathit{,}}{\mathit{\nu }}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\nu }}^{{\mathit{\beta }}_{\mathrm{d}}}\mathit{B}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathit{,\nu }\mathrm{\right)}}{\frac{\mathit{\partial B}\mathrm{\left(}{\mathit{T}}_{\mathrm{CMB}}\mathit{,\nu }\mathrm{\right)}}{\mathit{\partial T}}}\hspace{0.17em}\mathit{/}\hspace{0.17em}\frac{{\mathit{\nu }}^{{\mathit{\beta }}_{\mathrm{d}}}\mathit{B}\mathrm{\left(}{\mathit{T}}_{\mathrm{d}}\mathit{,}{\mathit{\nu }}_{\mathrm{0}}\mathrm{\right)}}{\frac{\mathit{\partial B}\mathrm{\left(}{\mathit{T}}_{\mathrm{CMB}}\mathit{,}{\mathit{\nu }}_{\mathrm{0}}\mathrm{\right)}}{\mathit{\partial T}}}\mathit{,}$(11)where the dust amplitude, A^Dust, is measured in units of μK², ν_i is the reference frequency for map i, ν₀ is a reference frequency, which is taken to be 143 GHz, B(T,ν) is the emission law of a black body with temperature T, and the dust colour-correction terms, ${\mathit{g}}_{\mathit{i}}^{\mathrm{Dust}}$, are computed by integrating the dust spectrum within the spectral band of each detector (set). We fix the frequency and angular scaling parameters to γ_d = 2.6, β_d = 1.6 and T_d = 18 K.

3.3. Poisson power from unresolved point sources

Unresolved galaxies contribute both shot noise and clustered power to the Planck maps. The Poisson contribution leads to a scale independent tem, C_ℓ = const. We model this power with a single amplitude parameter for each auto-spectrum (${\mathit{A}}_{\mathrm{100}}^{\mathrm{PS}}$, ${\mathit{A}}_{\mathrm{143}}^{\mathrm{PS}}$, and ${\mathit{A}}_{\mathrm{217}}^{\mathrm{PS}}$) and a cross correlation coefficient for each cross spectrum (${\mathit{A}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{PS}}\mathrm{=}{\mathit{r}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{{PS}}\sqrt{{\mathit{A}}_{\mathrm{143}}^{{PS}}{\mathit{A}}_{\mathrm{217}}^{{PS}}}$). These quantities are not of primary interest for cosmological results, so to avoid modelling error we do not separate the power into that sourced by “dusty" or “radio" galaxies (i.e., with increasing or decreasing brightness with frequency, respectively) as is done in the analysis of the ACT and SPT power spectra (Dunkley et al. 2011, 2013; Reichardt et al. 2012). We also make no assumptions about their coherence between frequencies.

The Poisson power is related to the flux density dN/dS via ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\mathrm{\int }\mathrm{d}\mathit{n̂}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{S}}_{\mathrm{cut}}\mathrm{\left(}\mathit{n̂}\mathrm{\right)}}\mathrm{d}\mathit{S}\hspace{0.17em}{\mathit{S}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathit{,}$(12)where we have explicitly introduced the Planck flux cut . Since Planck utilises a constant signal-to-noise cut, and the Planck noise varies significantly across the sky, this flux cut has a spatial dependence. Although this does not alter the shape of the Poisson term, extra care must be taken when comparing results ⁴ with models of dN/dS. In Sect. 7.3, we explore the consistency between the Poisson power recovered from the Planck power spectrum analysis and predictions from source count measurements.

3.4. Clustered power from unresolved point sources

Unresolved galaxies also contribute power because they trace large-scale structures. The mean flux from the radio galaxies is much smaller than that from the dusty galaxies, so only the dusty galaxies contribute a significant clustering term (Millea et al. 2012). The CIB clustering has been studied extensively, starting with Bond et al. (1986, 1991). Further theoretical investigation (Scott & White 1999; Haiman & Knox 2000) was stimulated by the detection of the infrared background in the COBE data (Puget et al. 1996; Fixsen et al. 1998), and the detection of bright “sub-millimetre” galaxies in SCUBA data (Hughes et al. 1998). Subsequently, the clustering has been detected at 160 microns (Lagache et al. 2007), at 250, 350 and 500 microns by the Balloon-borne Large Aperture Submillimeter Telescope (BLAST, Viero et al. 2009; Hajian et al. 2012) and at 217 GHz by SPT and ACT (Hall et al. 2010; Dunkley et al. 2011). Recent Planck measurements of the CIB (Planck Collaboration XVIII 2011) have extended the measurements at 217 GHz, 353 GHz, and 545 GHz to larger scales, and recent Herschel measurements (Viero et al. 2013) have improved on the BLAST measurements and extended them to smaller angular scales.

Rather than attempt to establish a physical model of the CIB, we adopt in this analysis a phenomenological model that captures the CIB uncertainties for both Planck and high-ℓ experiments. Our baseline model for the clustered CIB component is a power-law spectrum with a free spectral index, ${\mathrm{𝒟}}_{\mathit{\ell }}^{\mathrm{CIB}}\mathrm{\propto }{\mathit{\ell }}^{{\mathit{\gamma }}_{\mathrm{CIB}}^{}}\mathit{,}$ with an amplitude at each frequency, ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}$ and ${\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}$, and a cross-correlation between frequencies, ${\mathit{A}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{CIB}}\mathrm{=}{\mathit{r}}_{\mathrm{CIB}}\sqrt{\mathrm{(}}{\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}{\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}\mathrm{)}$. We assume that the CIB clustering power at 100 GHz is negligible.

3.5. Unresolved Sunyaev-Zeldovich effects

Based on analysis of ACT and SPT data, the thermal Sunyaev-Zeldovich (tSZ) contribution is expected to contribute approximately ${\mathrm{𝒟}}_{\mathit{\ell }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{3000}}^{\mathrm{tSZ}}\mathrm{~}\mathrm{9}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ at 100 GHz and ${\mathrm{𝒟}}_{\mathit{\ell }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{3000}}^{\mathrm{tSZ}}\mathrm{~}\mathrm{4}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ at 143 GHz (Reichardt et al. 2012; Dunkley et al. 2013; Sievers et al. 2013). The kinetic Sunyaev-Zeldovich (kSZ) effect is expected to have a similar, or smaller, contribution, with ${\mathrm{𝒟}}_{\mathit{\ell }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{3000}}^{\mathrm{kSZ}}\lesssim \mathrm{7}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$. In addition, theoretical arguments (Reichardt et al. 2012; Addison et al. 2012b) suggest that there should be a tSZ × CIB correlation that should contribute about the same order of magnitude as the kSZ term at 143 GHz.

For Planck, all of these SZ contributions are small in comparison to other unresolved foregrounds and are therefore poorly constrained by Planck data alone. Nevertheless, to eliminate biases in cosmological parameters (Millea et al. 2012; Zahn et al. 2005), we model their contributions, with appropriate constraints from higher resolution CMB experiments, using three templates.

First, for the thermal SZ effect we adopt the family of templates described by Efstathiou & Migliaccio (2012). These are based on the Komatsu & Seljak (2002) model, but use the “universal” X-ray electron pressure profile, P_e, of Arnaud et al. (2010) extrapolated to high redshift via ${\mathit{P}}_{\mathrm{e}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{\propto }\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{\Lambda }}\mathrm{\left)}\mathrm{\right(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{\Lambda }}{\mathrm{]}}^{\mathrm{4}\mathit{/}\mathrm{3}\mathrm{-}\mathit{ϵ}\mathit{/}\mathrm{2}}\mathit{.}$(13)Here, ϵ describes departures from self-similar evolution, and a value of ϵ = 0.5, which is adopted as the default for parameter estimation purposes, provides a good match to the results from recent hydrodynamical numerical simulations incorporating feedback processes (Battaglia et al. 2010, 2012). Figure 6 shows the tSZ templates for three values of ϵ; the template shape is not particularly sensitive to ϵ. We treat the (dimensionless) normalization of the tSZ template at 143 GHz as an adjustable parameter, ${\mathrm{𝒟}}_{\mathit{\ell }}^{\mathrm{tSZ}}\mathrm{=}{\mathit{A}}_{\mathrm{143}}^{\mathrm{tSZ}}{\mathrm{𝒟}}_{\mathit{\ell }}^{\mathrm{tSZ}\mathit{,}\mathrm{template}}\mathit{.}$(14)This parameter fixes the amplitude at 100 GHz via the frequency dependence of the tSZ effect, ${\mathit{C}}^{\mathrm{tSZ}}\mathrm{\propto }{\left(\mathit{x}\frac{{\mathrm{e}}^{\mathit{x}}\mathrm{+}\mathrm{1}}{{\mathrm{e}}^{\mathit{x}}\mathrm{-}\mathrm{1}}\mathrm{-}\mathrm{4}\right)}^{\mathrm{2}}\mathit{,} \mathit{x}\mathrm{=}\frac{\mathit{h\nu }}{\mathit{kT}}\mathrm{·}$(15)We neglect the tSZ at 217 GHz.

Second, for the kinetic SZ effect we adopt the template described by Trac et al. (2011), and as in Eq. (14) we treat the dimensionless amplitude of the template, A^kSZ, as a free parameter, $\begin{array}{ccc}{\mathrm{𝒟}}_{\mathit{\ell }}^{\mathrm{kSZ}}\mathrm{=}{\mathit{A}}^{\mathrm{kSZ}}{\mathrm{𝒟}}_{\mathit{\ell }}^{\mathrm{kSZ}\mathit{,}\mathrm{template}}\mathit{.}& & \end{array}$(16)Third and finally, for the cross-correlation between the thermal SZ component and the CIB we adopt the template described by Addison et al. (2012a). In this case, the amplitude is parametrised in terms of a single correlation coefficient, $\begin{array}{c}\\ \\ \end{array}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}$(17)These templates are plotted in Fig. 6, normalized to σ₈ = 0.8 and with ξ = 1.0 using a fiducial CIB amplitude. Note that with these parameters, the tSZ × CIB cross-spectrum approximately cancels the kSZ spectrum at 143 GHz.

As seen in Fig. 6, the SZ contributions are at the level of a few μK², which, although small, must be taken into account to assess inter-frequency residuals. However, one can see that these templates have similar shapes at multipoles ≲2000, and therefore they cannot be disentangled using Planck data alone. On the other hand, higher resolution experiments can break this degeneracy, and as shown in Planck Collaboration XVI (2014), the combination of Planck, ACT, and SPT, better constrains the amplitude of the thermal SZ effect. The ACT and SPT data at 150 GHz can be fitted to high, sub-μK², accuracy without kSZ and tSZ x CIB templates, yet we expect a kSZ contribution of at least the amplitude shown in Fig. 6, and larger if we account for patchy reionisation (see, e.g., Knox 2003, and references therein) and references therein). This implies a cancellation of the kSZ and tSZ x CIB contributions at 150 GHz (Addison et al. 2012b), as discussed in greater detail in Planck Collaboration XVI (2014).

4. Consistency checks and combined cross-spectra

The large number of cross spectra in a detector-by-detector power spectrum analysis allows for a number of internal consistency checks of the data. Within a frequency band, we expect to see exactly the same sky signals (primordial CMB, Galactic, and extra-galactic foregrounds), and so any intra-frequency residuals reflects instrumental systematics, for example beam errors, “gain” fluctuations, and band-pass mismatch. In contrast, inter-frequency residuals are harder to analyse because the sky signals vary with frequency. An accurate model of the unresolved foregrounds is therefore required to assess inter-frequency residuals. Furthermore, as we show below, the scatter caused by chance CMB–foreground cross-correlations can dominate the inter-frequency residuals. For a precision experiment such as Planck, where the power spectra are expected to be signal dominated over a wide multipole range, intra- and inter-frequency residuals provide a powerful way of assessing possible systematic errors. It is essential that contributions of systematic errors to both types of residual are small enough that they have negligible impact on cosmological parameter analysis.

Figure 7 shows a selection of temperature cross-spectra and estimates of the analytic covariance matrices, together with the best-fit cosmological model described in Sect. 5. Unresolved foregrounds have been subtracted using the best-fit foreground parameters of the model described in Sect. 3. The scatter varies substantially between cross-spectra, reflecting differences in the instrument noise and effective resolution of different detector combinations. The analytic error model summarized in Appendix A.2 is indicated, modified by the non-white noise correction. This model provides an excellent description of the scatter seen in the data, over the full multipole range shown in plots, with an accuracy of a few percent or better.

4.1. Intra-frequency residuals

In this section we analyse the intra-frequency residuals at 143 and 217 GHz. There are N_spec = 10 cross-spectra at 143 GHz and 15 cross-spectra at 217 GHz⁵. At each frequency, we solve for multiplicative (“effective” calibration) coefficients, y_i, that minimise ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{\ell }}\sum _{\mathit{ij,j}\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathit{i}}\mathrm{(}{\mathit{y}}_{\mathit{i}}{\mathit{y}}_{\mathit{j}}\mathit{Ĉ}\begin{array}{c}\mathit{ij}\\ \mathit{\ell }\end{array}\mathrm{-}\mathrm{⟨}\mathit{Ĉ}\mathit{\ell }\mathrm{⟩}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(18)where $\mathrm{⟨}\mathit{Ĉ}\mathit{\ell }\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{spec}}}\sum _{\mathit{ij,j}\hspace{0.17em}\mathit{>}\hspace{0.17em}\mathit{i}}{\mathit{y}}_{\mathit{i}}{\mathit{y}}_{\mathit{j}}\mathit{Ĉ}\begin{array}{c}\mathit{ij}\\ \mathit{\ell }\end{array}\mathit{,}$(19)subject to the constraint that y₁ = 1 (where i = 1 corresponds to detector 5 at 143 GHz and detector 1 at 217 GHz). Note that the power spectra in Eq. (18) and (19) are corrected for beam transfer functions. To minimise the possible impact of beam errors and noise, we restrict the sum in Eq. (18) to the multipole range 50 ≤ ℓ ≤ 500 where the spectra are signal dominated. Numerical values for the calibration coefficients are given in Table 3, using mask CL31. The calibration factors are insensitive to the choice of mask or multipole range.

The results of Table 3 show that effective calibration factors of ~0.2% are quite typical for HFI maps, in the 100 − 217 GHz frequency range. These recalibrations are of the order of magnitude of the statistical errors of the calibrations on dipole (see Planck Collaboration VIII 2014, Table 2). Note that the data are corrected for individual bolometer time transfer functions (TTFs; Planck Collaboration VII 2014). For each detector, the TTF model is tuned to minimise survey differences and by construction normalized to unity at the spin frequency of the satellite (0.01666 Hz) to preserve the dipole calibration. The consistency of intra-frequency power spectrum residuals therefore provides a test of the consistency of the TTFs in addition to the beam transfer functions.

Figure 8 shows the remarkable consistency of the power spectra at each frequency. The upper panels show the spectra corrected for the beam and effective calibration, together with the mean cross spectra. The lower panels show the dispersion around the mean. Up to ℓ ≃ 700 the Δℓ = 31 binned cross-spectra show a dispersion in their values typically less than expected from ~0.2% calibration variations. This excess scatter (compared to the noise-induced expected scatter) has negligible impact on cosmological parameter analysis.

The residuals of the cross spectra in band averages of width Δℓ ~ 61 are shown in Fig. 9, before and after correction for the effective intra-frequency calibrations. The reduction in scatter after correction is evident at ℓ ≲ 500, and the residual scatter is consistent with instrument noise and beam errors. At 217 GHz, beam errors dominate over noise at multipoles ≲1000. There is no evidence that the excess scatter is caused by a small number of “anomalous” detectors.

4.2. Inter-frequency residuals

The results of the previous section show that the intra-frequency cross-spectra between detector/detector sets are consistent to within a few μK² at multipoles ℓ ≲ 1000. In a likelihood analysis, there is therefore little loss of information in compressing the power spectra for each distinct frequency combination, as opposed to retaining the spectra for each map pair. This compression greatly reduces the size of the data vector and its covariance matrix, and speeds up the likelihood computation at high multipoles. In this section we inter-compare the residuals of these compressed power spectra.

One might naïvely expect that with accurate foreground modelling, the inter-frequency residuals in the signal dominated regime should be reduced to levels comparable to those seen in the intra-frequency comparisons described in the previous section. This is incorrect. Figure 10 shows power spectrum differences between the cosmologically significant spectra for Planck at high multipoles (143 × 143, 143 × 217, 217 × 217). In this figure, which is independent of the cosmological model, the best-fit unresolved foreground model has been subtracted from each spectrum, and relative calibration factors have been applied. Residual beam, calibration and unresolved foreground errors would show up in this figure as large-scale smooth residuals.

In fact, we see small-scale residuals at multipoles ℓ ≲ 800, which are larger than expected from instrument noise alone. For example, the 143 × 143 − 217 × 217 spectrum in the upper panel of Fig. 10 has a χ² of 35.6 for the 21 points in the multipole range 100 ≤ ℓ ≤ 800, an approximately 2.3σ excess. (But note that this effect would become more significant on the full mission data.) This excess scatter arises from the chance CMB-foreground cross-correlations. Even if the foreground contamination is much smaller than the CMB, chance cross-correlations can produce scatter in the inter-frequency power spectra that dominates over instrument noise in the signal dominated regime. We develop a model for these CMB-foreground cross-correlations in Appendix C and show it can explain quantitatively the observed inter-frequency scatter, and the scaling of this scatter with frequency, at multipoles ℓ ≲ 500.

At high enough multipoles, instrument noise, beam errors, and errors in foreground modelling dominate the inter-frequency residuals. A complete analysis of inter-frequency residuals therefore requires the full likelihood machinery and MCMC analysis to determine foreground, beam and calibration parameters. We therefore revisit the inter-frequency residuals in the following sections.

5. Reference results of the high-ℓ likelihood

In this section we study the high-ℓCamSpec likelihood, and present the power spectrum and parameters derived from this baseline likelihood for the basic six-parameter ΛCDM model. In order to break the well-known degeneracy between the optical depth, τ, and the scalar index of scalar perturbations, n_s, we adopt a Gaussian prior on τ (inspired from WMAP7 data, i.e., 0.088 ± 0.015, see Komatsu et al. 2011) instead of the low-ℓ likelihood at ℓ < 50. We return to the global Planck results after introducing the low-ℓ likelihood.

We choose separate masks for each frequency map to minimise Galactic foreground emission. First, since the HFI data are signal-dominated at ℓ ≲ 500, and diffuse Galactic emission is low at 100 GHz outside the CL49 mask, there is little to be gained from analysing the dustier 217 GHz with the same region; it contains no new information about the primordial CMB. At higher multipoles, ℓ ≳ 500, we use the CL31 mask, optimizing the sky coverage while ensuring a low amplitude of small-scale Galactic emission relative to the isotropic unresolved foregrounds and the primordial CMB. In addition, we tune the multipole range for each frequency to mitigate Galactic foreground contamination and beam errors.

The choices of masks and angular ranges used in the high-ℓ likelihood are summarized in Table 4, together with basic χ² statistics with respect to the minimal ΛCDM model per cross-spectrum and combined. The 100 GHz cross-spectrum is computed over the largest sky fraction, a total of 49% of the sky, and measures the largest scales. On the other hand, it has lowest resolution, and it is therefore only used for ℓ ≤ 1200. The 143 GHz cross-spectrum has higher resolution, and is used for ℓ ≤ 2000. Finally, the 217 GHz cross-spectrum has the highest resolution, but also the most Galactic dust contamination, and is therefore evaluated from only 31% of the sky, but including an angular range of 500 ≤ ℓ ≤ 2500.

Given these masks and angular ranges, we compute the angular power spectra and covariance matrices, and construct the CamSpec likelihood. The angular power spectra for each frequency combination are shown in Fig. 11, and compared to spectra derived from the 70 GHz and 353 GHz Planck maps.

We use the likelihood to estimate six ΛCDM cosmological parameters, together with a set of 14 nuisance parameters (11 foreground parameters described in Sect. 3, two relative calibration parameters, and one beam error parameter⁶, see Eq. (5) and Appendix A.7 for details). Tables 5 and 6 summarize these parameters and the associated priors⁷. Apart from the beam eigenmode amplitude and calibration factors, we adopt uniform priors. To map out the corresponding posterior distributions we use the methods described in Planck Collaboration XVI (2014), and the resulting marginal distributions are shown in Fig. 12. Note that on the parameters A^tSZ, A^kSZ and ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}$ we are using larger prior ranges as compared to Planck Collaboration XVI (2014).

Figure 12 shows the strong constraining power of the Planck data, but also highlights some of the limitations of a “Planck -alone” analysis. The thermal SZ amplitude provides a good example; the distribution is broad, and the “best fit” value is excluded by the ACT and SPT high resolution CMB experiments (Reichardt et al. 2012). For the CIB amplitudes, the upper bound on e.g., ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}$ is significantly weaker than the ACT and SPT constraints. To accurately estimate the foreground parameters at the ≲μK² level, we need to supplement the Planck power spectra with temperature data from ACT and SPT, as described in Planck Collaboration XVI (2014). The fiducial model and foreground parameters used in the CamSpec likelihood are therefore derived from a joint Planck+ACT+SPT analysis and is not based on the parameters listed in Table 8. In the rest of this section, we use the parameters of Table 8, which is derived using only Planck data, to discuss inter-frequency residuals.

Figure 13 shows the foreground residuals and total residuals after removing the best-fit foreground model for all spectra (including the 100 × 143 and 100 × 217 spectra, which are not used in the CamSpec likelihood). The first point to note here is that Planck has a limited ability to disentangle foregrounds. While the Planck data constrain the Poisson point source amplitudes at each frequency, as well as the CIB amplitude at 217 GHz (which dominates over the Poisson point source amplitude over much of the multipole range), they have only marginal sensitivity to the tSZ amplitude in the 100 × 100 spectrum, though the thermal SZ is strongly degenerate with the Poisson point source amplitude. The remaining foreground parameters are highly degenerate. For Planck alone, these minor foreground contributions combine to absorb inter-frequency residuals.

Pairs of spectra are compared in Fig. 14, averaged over bands of width Δℓ = 31 below ℓ ≲ 2000 and wider bands above 2000. The error bars show the diagonals of the covariance matrices of these averages, but it is important to note that these points are correlated (neighbouring points are actually anti-correlated at the few percent level) even with bin widths as large as these because of the sky masks. This comparison shows that each of the spectra used in the CamSpec likelihood is consistent with the best-fit theoretical spectrum to high accuracy. In fact, each spectrum can be used to form a likelihood, and each gives a reduced χ² close to unity (see Table 4). Thus, the six parameter ΛCDM model provides an excellent fit to the Planck high-ℓ power spectra at all frequencies between 100 and 217 GHz.

Figure 15 shows our maximum likelihood primary CMB spectrum, together with the best-fit theoretical spectrum. The residuals with respect to the model are shown in the lower panel. The error bars are computed from the diagonal components of the band-averaged covariance matrix. The binning scheme is the same as in Fig. 14.

Finally, in Fig. 16 we zoom in on this spectrum in four multipole ranges using finer binning. The correlated fluctuations seen in this figure are mask-induced, and perfectly compatible with the six parameter ΛCDM model. Features such as the “bite” missing from the third peak at ℓ ~ 800 and the oscillatory features in the range 1300 ≲ ℓ ≲ 1500 are in excellent agreement with what we expect from our covariance matrices and from simulations; see Appendix A.4 for a few specific examples.

6. Accuracy assessment of the high-ℓ likelihoods

In this section we compare the power spectra and likelihoods derived using our two independent methods, and test these likelihoods using full Planck simulations.

6.1. Comparison of the Plik and CamSpec likelihoods

To allow a more direct comparison between the CamSpec and Plik likelihoods, we use the same frequency cross-spectra for both codes in the following, i.e., we discard the 100 × 143 and 100 × 217  GHz frequency combinations from the default setting of the Plik likelihood. To achieve this, we modify the Plik likelihood to use the fiducial Gaussian approximation instead of the Kullback divergence. On the other hand, while we use the same multipole coverage, 100 ≤ ℓ ≤ 2500, and only one mask (CL39) for all cross-spectra for Plik, we still use multipole ranges and masks as defined in Table 4 for CamSpec. In addition, we perform one Plik analysis with the CL49 mask, which matches the CamSpec mask at 100 GHz.

The left column of Fig. 17 shows the differences between the Plik power spectra, adopting the above validation settings, and the corresponding best-fit model. The right column shows the total spectra decomposed into cosmological and foreground components. The residuals are similar to those in Fig. 13 in the sense that they do not show any evidence of biases, except for some excess power in the 217 × 217 GHz spectra at small scales, where foreground modelling has the highest impact. At scales ℓ ≲ 1500 (ℓ ≲ 1000 for the 100 × 100 GHz cross spectra), the residuals are coherent between cross spectra as they are computed with the same Galactic mask, and the residuals are dominated by (CMB-induced) cosmic variance. At smaller scales (ℓ ≳ 1500) the residuals are dominated by noise and become uncorrelated.

In Fig. 18 we show the CMB power spectrum recovered by Plik  estimated by removing the best-fit foreground amplitudes from each cross-spectrum and computing their optimally weighted average, and the corresponding difference with respect to the best-fit ΛCDM model. The large scatter at low multipoles is expected, due to cosmic variance. The residual scatter at higher multipoles is at the ± 20 μK² level, demonstrating the good fit provided by the sum of the ΛCDM and foreground models. These CMB residuals can be compared to the CamSpec inverse-covariance weighted CMB residuals shown in Fig. 15, which are of the same order of magnitude. Thus, the Planck likelihood fit to the ΛCDM model is robust with respect to the detailed shape of the likelihood, as quantified in terms of power spectrum residuals.

The ΛCDM parameter constraints derived from the two likelihoods are shown in Fig. 19, while Fig. 20 shows the foreground parameters. For the cosmological parameters, the agreement between the two likelihoods is excellent: when the CL49 mask is adopted for the Plik likelihood, which is also used by CamSpec at 100 GHz, all parameters agree to 0.2σ in terms of maximum posterior values. We also see that the widths of the distributions are quite similar, with the Plik ones slightly broader than the CamSpec ones. Significantly larger differences are seen for the foreground parameters.

These differences can be understood as follows: we use the CamSpec likelihood with Galactic mask CL49 for the 100 × 100 GHz spectra, to minimise the cosmic variance in the low to intermediate multipole range (ℓ ≲ 1200), taking advantage of the low level of Galactic emission in this channel. At higher multipoles we use the more conservative Galactic mask CL31 for both the 143 and 217 GHz channels, at the price of a higher variance. However, in the specific case of the ΛCDM model considered here, most of the constraints on cosmological parameters come from relatively modest multipoles, ℓ < 1500, rather than from the damping tail. This explains why when we repeat the Plik analysis, enlarging the sky area from Galactic mask CL39 to CL49, we find parameter distributions in good agreement with those of CamSpec.

Figure 20 shows the foreground parameters estimated from both Plik and CamSpec. We consider the case for Planck data alone, and with the inclusion of data from ACT and SPT. We also impose, for CamSpec, a Gaussian prior of 0.7 ± 0.2 on the CIB slope parameter, γ^CIB. We find that the upper bounds on the CIB amplitude at 143 GHz, ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}$, and on the SZ amplitudes (both thermal and kinetic, A^tSZ and A^kSZ) are in good agreement using Planck data alone, but we see differences (~1.5σ) in the CIB and Poisson amplitudes at 217 GHz (${\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}$, ${\mathrm{A}}_{\mathrm{217}}^{\mathrm{PS}}$) as well as a difference in the CIB correlation coefficient, ${\mathit{r}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{CIB}}$.

We understand this effect in the following way: CamSpec uses a more limited multipole range and a more conservative mask at 217 GHz than Plik, in order to minimise the Galactic emission in this channel. This enhances the degeneracy between the foreground parameters at 217 GHz, and enhances the sensitivity to possible deviations of the CIB power spectrum from the pure power law assumed here. This is artificially enhanced by normalizing the components at ℓ = 3000, which is more suitable for high resolution experiments than for Planck. Let us stress that, despite the disagreement of the precise decomposition into the different physical components, the sum of the foreground contributions at 217 GHz in Plik and CamSpec is in good agreement. And when ACT and SPT data are added, these differences are largely reduced, as shown in Fig. 20.

Are these differences important for cosmology? To address this question, we examine the covariance between cosmological and foreground parameters. Figure 21 shows the correlation matrix for all the estimated parameters. The basic ΛCDM parameters have well-known correlations: the scalar spectral index, n_s, is anti-correlated with both the amplitude A_s and the dark matter density Ω_ch², which is itself correlated with the amplitude. Within the foreground parameters, there are strong correlations between the Galactic dust amplitude, the CIB amplitudes, and the SZ amplitude, as well as between point source amplitudes at different frequencies. These correlations result from the conservative foreground model adopted here, where all amplitudes of the CIB and Poisson contributions are left free to vary in each frequency pair; this choice results in small residuals in the fits, at the price of partial degeneracies between the foreground parameters.

Despite this conservative foreground model, there is a small correlation between cosmological and foreground parameters, the strongest effect being a 34% anti-correlation between the kinetic SZ amplitude and the scalar spectral index. The kinetic SZ power spectrum amplitude is positive, so marginalising over it affects the peak value of n_s despite the fact that A^kSZ is not significantly different from zero. The addition of smaller scale data from ACT and SPT helps to break this degeneracy. In addition, as can be seen in Fig. 20, the posterior distributions of CamSpec and Plik for A^kSZ are in good agreement, showing the stability of the cosmological parameters to the likelihood method.

In Sect. 7, we further explore the stability of the cosmological and foreground distributions to technical choices made in the likelihood and data selections.

6.2. Comparison to simulations

We now investigate the precision and accuracy with which cosmological parameters, and foreground parameters to a lesser extent, can be recovered from realistically simulated data. Here we compare the posterior distributions of cosmological and foreground parameters, together with calibration and beam error parameters, inferred using the Plik likelihood, with the input values of a set of simulations, referred to as “Full Focal Plane” (FFP6). The signals in these simulations are based on the “Planck Sky Model” (Delabrouille et al. 2013), which includes a detailed model of the astrophysical emission, both Galactic and extragalactic, at the Planck frequencies. The simulations also reproduce in detail the main instrumental systematic effects of Planck, including correlated time-line noise, instrumental pointing, flags, anisotropic detector beams, and spectral bandpasses. One thousand CMB and noise realisations were generated using the same foreground emission, and a hundred realisations were performed at the level of different detector sets. These simulations are described further in Planck Collaboration (2013).

In order to test for any methodological bias, we estimate cosmological parameters from 100 FFP6 simulations using the Plik likelihood. They consist of random realisations of CMB anisotropies and noise, superimposed with a single realisation of a frequency-dependent foreground template. We assume the foreground power spectra to be known exactly, and include their additional power as a constant component in our model. We do not estimate the foreground parameters in this case. The goal of this test is to demonstrate the reliability of our analysis pipeline and explore the effects of the noise and CMB-foreground chance correlations on parameter estimation.

To compare our results to the simulation inputs, we also estimate cosmological parameters from the 100 CMB realisations, without including noise or foregrounds. In the likelihood evaluations we down-weight the high-ℓ part of the power spectrum as if there were noise at the level of the noise simulations. A direct comparison of the derived parameters from this CMB-only analysis and the full simulation (that contains noise and foreground emission) allows us to remove the scatter introduced by cosmic variance.

The result of the analysis is shown in Fig. 22, where we plot the distribution of the difference of the mean estimated parameters, between the CMB-only simulation and the noisy CMB simulation with foregrounds, in units of the standard deviation of each individual distribution. An unbiased pipeline will give a difference consistent with zero. Averaged over 100 simulations, only the bias on Ω_ch² is statistically significant, but is <0.3σ. An interesting point to note from these histograms is that the distributions are rather wide, mostly around 0.5σ, while by construction, we only studied the effect of noise and CMB-foreground chance correlations in this test.

The FFP6 simulations can also be used to assess the sensitivity to foreground modelling errors. In Fig. 23, we compare the posterior marginal distributions of cosmological parameters, estimated with different assumptions about the foreground model. The blue lines, which correspond to the analysis of CMB-only simulations (but accounting for the noise covariance in the likelihood) correspond to the idealised case where the foregrounds play no role, and where the noise-induced variance has been averaged. The red, thick lines show results obtained when marginalising over the parameters of the model of foregrounds that is applied to the Planck data (see Sect. 3), with a fixed value of the CIB spectral index γ^CIB. Purple and grey lines show respectively the effect of leaving γ^CIB free when marginalising, or fixing it to the displaced value of 0.4, more than 2σ away from the peak posterior. The green lines show the effect of leaving γ^Dust (spectral index of the Galactic dust emission) free in the marginalisation.

The distributions are all in reasonable agreement with the input parameters of the simulation. In addition, we see that varying assumptions on the parameters of the foreground model (red, green, purple, and yellow lines) have negligible impact on the recovered cosmological parameters. Finally, the broadening of the posteriors between the CMB-only exploration and the full case (including noise random realisations and foregrounds) is expected, as the latter includes all the sources of variance, including CMB-foreground and CMB-noise chance correlations.

It is worth noting that the FFP6 foreground simulations, based on extrapolations of existing observations, cannot be described by the simple foreground model used in the likelihood analyses. The upper bounds on biases introduced by a possible mismatch between the simulated foreground templates and the model used in the analysis, inferred from Fig. 22, should be representative of the Planck data analysis. The negligible impact of the various assumptions made on the foreground model parameters of Fig. 23 confirms that the cosmological parameter estimations are robust to details of the foreground model.

7. Consistency checks

In this section we investigate the stability of the distributions of cosmological and foreground parameters. The technical choices made in constructing the high-ℓ likelihood fall into three broad categories. The first category covers internal parameter choices that leave the data selection unchanged. This includes choices such as the binning strategy, marginalising or not over calibration and beam errors, and the description of the noise model. The second category includes variations in the data selection, such as the multipole range used, and the choice of masks, or of spectra. The final category accounts for variation in the foreground model. We perform a suite of tests to investigate the impact of these choices on parameters. We use the Plik likelihood, and all tests are compared to the baseline Plik spectra. Most of the results can be summarized by “whisker plots”, which compare the main properties of the posterior distribution of the cosmological and foreground parameters. More detailed results are reported in Appendix D.1.

In this section we also compare our estimated cosmological parameters to those derived from spectra computed from the LFI 70 GHz channel. We additionally check the consistency of parameters with results obtained using the power spectrum of CMB maps derived by component separation methods (described in Planck Collaboration XII 2014) that use Planck data at all frequencies. This battery of tests demonstrates the stability of the inferred cosmological parameters.

A final test is to compare the predicted polarisation spectrum of the best fitting ΛCDM model with spectra measured from Planck. As discussed in Planck Collaboration I (2014) and Planck Collaboration VI (2014), the Planck polarisation data is not yet used in our cosmological analysis, as further tests must be performed, but the current results increase our confidence in the robustness of the high-ℓ temperature likelihood.

7.1. Impact of technical choices for fixed data selection

Here we consider three changes: (1) fixing the inter-frequency calibration and beam errors to the best-fit values, rather than marginalising; (2) including sub-pixel effects; and (3) including noise correlation between detectors. Figure 24 shows the corresponding impact on parameters.

The effect of fixing the calibration and beam errors is negligible on most cosmological parameters within the six parameter ΛCDM model, with the exception of n_s where we see a 0.16σ shift. There is a bigger effect on the foreground parameters due to their partial degeneracy, and their sub-dominant contribution to the total power. There is only a small correlation between the cosmological parameters and the calibration coefficients, so marginalising or fixing their value has little impact on the cosmology.

The calibration coefficients are, however, strongly correlated with each other, in particular at 217 GHz, since they are also significantly correlated with e.g., the CIB amplitudes in the 217 × 217 and 217 × 143 GHz spectra. This is important to keep in mind when comparing the calibration estimates obtained here with those obtained from the CMB dipole in the HFI data processing paper (Planck Collaboration VI 2014). Nevertheless, as can be seen in Fig. 26, the peak posterior values of the relative calibration coefficients are found to differ from 1 at most at the few parts per thousand for all the 13 detectors sets involved, in agreement with the estimates of the calibration accuracy of the maps (Planck Collaboration VI 2014), although a wide flat prior has been applied on these coefficients. The same test applied on simulations with no beam or calibration errors shows how well this test is passed. This confirms that the deviations found at the 0.1% are significantly detected, and it is important to show that these deviations have little impact on the cosmology.

The estimated values of the beam errors do not imply that extra beam corrections are required. A comparison of the prior and posterior distributions suggests that we have quite conservatively estimated the uncertainties from the beam determination. We keep this conservative approach in CamSpec, in which we marginalise analytically over all beam eigenmodes except for the dominant 100 GHz mode, ${\mathit{\beta }}_{\mathrm{1}}^{\mathrm{1}}$, which we sample directly.

As an extended test, we investigate the effect of possible errors in the beam transfer function when the helium abundance Y_P is also allowed to vary freely (i.e., without imposing constraints from Big Bang Nucleosynthesis), as it has a larger effect on the small scale spectrum. Varying this parameter leads to a substantial broadening of the posterior distributions for Ω_bh², θ_MC, and n_s, as can be seen in Fig. 27. We confirm that marginalising over calibration and beam errors has a small impact on all cosmological parameters, including Y_P. On the other hand, it has a somewhat larger impact on some of the foreground parameters (see Appendix D.1).

We then investigate the impact of two sub-dominant effects: the “sub-pixel effect” and the possible presence of a correlation in the noise between detectors or detector sets. The sub-pixel effect has a convolving effect on the power spectra that is similar to gravitational lensing of the CMB, but is purely a result of the Planck scanning strategy and the map-making procedure (Planck Collaboration VII 2014). The scanning strategy on rings with very low nutation levels results in the centroid of the samples being slightly shifted from the pixel centres; however, the map-making algorithm assigns the mean value of samples in the pixel to the position at the centre of the pixel. This has a non-diagonal effect on the power spectra, but the correction can be computed given the estimated power spectra for a given data selection, and recast into an additive, fixed component of the model covariance matrix.

The possible noise correlation between detectors may appear due to factors such as common residual thermal fluctuations, electronic chain noise, or cosmic ray showers. To build a model of this correlated component, we compute the cross-power spectra between detectors of difference maps that are free of signal. This procedure should capture all correlations on time scales shorter than half a ring’s observation. These estimates are noisy, so we compute an average amplitude of the correlation for ℓ ≥ 1000.

We find that the impact on cosmological parameters of both these effects is negligible, with less than a 0.1σ variation on the mean posterior values.

7.2. Impact of data selection

Here we consider three changes: (1) varying the angular range used in the likelihood; (2) varying the Galactic mask; and (3) discarding individual frequency channels. These are expected to result in changes in the parameter distributions due to the fact that we are changing the input data. The corresponding impact on parameter distributions are illustrated in Fig. 24.

We first vary the maximum and minimum multipole. Using ℓ_max = 1008 gives parameter distributions of similar width to those obtained by WMAP. We find that all basic ΛCDM cosmological parameters have converged by ℓ_max ≃ 1500, since no parameters are specifically sensitive to the damping tail. The convergence of the posteriors on foreground parameters is, as expected, slower when increasing ℓ_max, as most of them are dominant at small scales. The Galactic dust normalization decreases with ℓ_max due to its correlation with the CIB components. Similarly, changing ℓ_min from 100 to 50 has a negligible effect on cosmological parameters, and mostly affects the determination of the Galactic dust amplitude, which decreases for ℓ_min as it is better measured on large scales; its correlations with the CIB clustered and Poisson contributions explain the slight variations in the corresponding parameters ${\mathrm{A}}_{\mathrm{217}}^{\mathrm{CIB}}$, γ^CIB, and ${\mathrm{A}}_{\mathrm{217}}^{\mathrm{PS}}$.

We then investigate the impact of varying the Galactic mask, from the most conservative (CL20) to the least conservative (CL49), for fixed multipole range 100 ≤ ℓ ≤ 2508. As expected, the errors decrease as the sky fraction increases. From CL20 to CL39, cosmological and foreground parameters are stable. The foreground parameters change significantly, however, when we use CL49 for all channels, showing that our foreground model is unable to properly fit the data: the clearest sign of this failure is the unphysically low value of the CIB spectral index (γ^CIB), indicating that our CIB component determination is getting contaminated by a dust-like component with a steeper angular power spectrum than the CIB, but shallower than that of our (single) Galactic dust component. At low Galactic latitudes, the presence of compact Galactic sources leads to a flatter angular power spectrum than that of (high-latitude) diffuse thermal dust, and this likely affects the CIB determination. This justifies the choice made in CamSpecto use a conservative masking strategy (CL31) for the 143 and 217 GHz channels.

Next, we remove one frequency channel at a time. Results change by less than 0.5σexcept when removing the 217 GHz channel. Removing the 217 GHz channel retains only 21 cross-spectra out of the 78. As a finer test, we consider survey-difference tests, which were performed in the HFI data processing paper (Planck Collaboration VI 2014) at each frequency, for all combinations of the input maps used in the likelihood analysis (two at 100 GHz, five at 143 GHz, and six at 217 GHz), and for two survey differences: Survey 1 – Survey 2 and Survey 1 – Survey 3. It was found there that in addition to the 217-ds1 × 217-ds2 cross spectrum, only two other cross-spectra failed this test, namely 217-1 × 217-ds2 and 217-1 × 217-3. We therefore consider here a data selection where only the (other) “best 217 GHz” spectra are retained in the analysis. Figures 24 and 25 show that this has no discernible influence on the determination of all basic ΛCDM cosmological parameters (but for a slight broadening of the posterior distribution), while the main foreground parameter change is a 1σ shift of the point source amplitude at 217 GHz.

We also note that since the public release of the data in March 2013, we have found evidence that the “ℓ ~ 1800 dip” in some 217 GHz detector cross-spectra (see Fig. 17) is stronger in the first six-month survey than in subsequent surveys and that its amplitude may be reduced by additional data flagging targeted at EMI-EMC interference from the 4 K cooler drive and read-out electronics. This dip is therefore likely to be a (small) residual systematic effect in the data, but one, which has little impact on cosmological parameter determination in the standard model. Planck Collaboration XVI (2014) extends the verification to other cosmological models and data selections.

7.3. Testing the extragalactic sources model

In this section we describe various tests of the validity of our model for the extragalactic sources. Further tests are also reported in Appendix B of Planck Collaboration XVI (2014).

7.3.1. Poisson power from extragalactic sources

Here we check that the Poisson power estimated in the likelihood, which comes from sources below Planck’s detection threshold, is roughly consistent with the level expected given number counts of detected galaxies. Figure 28 shows source counts from Planck (Planck Collaboration Int. VII 2013) at 100, 143, and 217GHz derived from the Planck Early Release Compact Source Catalogue (Planck Collaboration VII 2011; Planck Collaboration XIII 2011). At 143 and 217 GHz we also show the source counts from SPT as reported in Vieira et al. (2010) at 150 and 220 GHz, and from ACT (Marriage et al. 2011) at 150 GHz. The models of de Zotti et al. (2005) and Tucci et al. (2011) are also shown and are discussed in Planck Collaboration Int. VII (2013).

Planck Collaboration Int. VII (2013) use spectral information to separate the sources into “synchrotron” and “dusty” sources, and show that the counts at 100–217 GHz are dominated by synchrotron sources at flux densities above ~400 mJy. Vieira et al. (2010) performed a similar separation. The counts at 150 GHz are dominated by synchrotron sources at flux densities S> 10 mJy, but dusty galaxies contribute roughly equally at 220 GHz at flux densities ≲30 mJy (Vieira et al. 2010; Hall et al. 2010). The ACT counts have not been separated according to spectral type, but should be dominated by radio sources at these flux densities.

Figure 28 show models fit to the counts using the function ${\mathit{S}}^{\mathrm{5}\mathit{/}\mathrm{2}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{=}\frac{\mathit{A}{\mathit{S}}^{\mathrm{5}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathit{S}\mathit{/}{\mathit{S}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\alpha }}\mathrm{+}\mathrm{(}\mathit{S}\mathit{/}{\mathit{S}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}}\mathrm{+}\mathit{B}\left(\mathrm{1}\mathrm{-}\exp \left(\frac{\mathit{S}}{{\mathit{S}}_{\mathrm{1}}}\right)\right)\mathit{,}$(20)where A, S₀, α, β, B, and S₁ are free parameters. The best-fit values of these parameters are given in Table 7.

Given this model, and given the approximate flux cut applied to the Planck maps, the expected contribution of radio sources to the Planck power spectra, at flux densities smaller than 400, 350, and 225 ± 50 mJy at 100, 143, and 217 GHz, are 8.47 ± 1, 6.05 ± 0.8, and 3.10 ± 0.7 Jy²/sr. The contribution of unresolved infrared galaxies to the power spectra is not negligible. They are expected to dominate at 217 GHz, even if they are sub-dominant in the Planck counts. Indeed, faint IR galaxies create a “bump” in the S^{5 / 2}dN/ dS distribution, below the detection limit of ACT or SPT. This bump is seen at higher frequencies, e.g., with the Herschel SPIRE instrument (see Planck Collaboration Int. VII 2013, for details).

This bump of infrared galaxies has not been measured at frequencies of 217 GHz and below. However, measurements with the AzTEC telescope at 1.1mm (270 GHz, Scott et al. 2012) can be used to extrapolate the counts down to 217 GHz. This leads to a predicted peak in the number counts (S^{5 / 2}dN/ dS) around 1.4 mJy at a level of 190 Jy^1.5/ sr, somewhat higher than the values in Hall et al. (2010). The corresponding contribution of infrared galaxies to the power spectrum is estimated in Planck Collaboration XVIII (2011) to be 16 Jy²/ sr at 217 GHz  but with significant uncertainty.

Summing the expected contributions from radio and IR galaxies, we estimate the following values for D₃₀₀₀ for Planck: 200, 75, and 120 μK² at 100, 143, and 217 GHz respectively. These predictions are much less certain at 217 GHz due to the absence of infrared galaxy counts at this frequency.

7.3.2. Clustered power in CIB fluctuations

Here we check the consistency of the estimated clustered power in CIB fluctuations. As already noted, from the CamSpec and Plik likelihoods we find only an upper limit on the clustered CIB power at 143 GHz. With CamSpec we detect clustered power at 217 GHz with ${\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}\mathrm{=}\mathrm{32}\mathrm{\pm }\mathrm{10}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$, and Poisson power with ${\mathit{A}}_{\mathrm{217}}^{\mathrm{PS}}\mathrm{=}\mathrm{92}\mathrm{\pm }\mathrm{22}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$. This Poisson power is dominated by the CIB fluctuations. For the Plik likelihood, we have ${\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}\mathrm{=}\mathrm{49}\mathrm{\pm }\mathrm{7}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ and ${\mathit{A}}_{\mathrm{217}}^{\mathrm{PS}}\mathrm{=}\mathrm{58}\mathrm{\pm }\mathrm{19}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$. The sum of the CIB power at 217 GHz, and at pivot scale ℓ = 3000, is in the range ≈105–125 μK²⁸.

We compare this level to the measurements by the ACT and SPT experiments, which probe higher angular resolution. Fitting a common model to the ACT power spectra from Das et al. (2014), and the SPT spectra from Keisler et al. (2011); Reichardt et al. (2012), the analysis in Dunkley et al. (2013) finds ${\mathit{A}}_{\mathrm{219.6}}^{\mathrm{CIB}}\mathrm{=}\mathrm{54}\mathrm{\pm }\mathrm{16}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ for the CIB clustered component (and ${\mathit{A}}_{\mathrm{219.6}}^{\mathrm{PS}\mathit{,}\mathrm{CIB}}\mathrm{=}\mathrm{78}\mathrm{\pm }\mathrm{12}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ for the CIB Poisson component) at effective frequency 219.6 GHz. For SPT, the clustered level is ${\mathit{A}}_{\mathrm{219.6}}^{\mathrm{CIB}}\mathrm{=}\mathrm{59}\mathrm{\pm }\mathrm{12}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$ (and Poisson ${\mathit{A}}_{\mathrm{219.6}}^{\mathrm{PS}\mathit{,}\mathrm{CIB}}\mathrm{=}\mathrm{69}\mathrm{\pm }\mathrm{10}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$), also at an effective frequency of 219.6 GHz. This is consistent with the SPT analysis in Reichardt et al. (2012). These are estimated assuming γ^CIB = 0.8 and ${\mathit{r}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{CIB}}\mathrm{=}\mathrm{1}$, and that the CIB emission can be modelled with frequency as a modified blackbody, following Addison et al. (2012a).

The total CIB signal seen by Planck, extrapolated to ℓ = 3000 scales, is therefore consistent with the ACT and SPT observations, but given the limited angular range of Planck, the clustered and Poisson part are degenerate. This motivates us to include the ACT and SPT data in many of our cosmological analyses.

When combining Planck, ACT and SPT data together, using the same foreground model (except for Poisson power, which depends on the respective flux cuts of the experiments), both CamSpec and Plik give ${\mathit{A}}_{\mathrm{217}}^{\mathrm{CIB}}\mathrm{=}\mathrm{50}\mathrm{\pm }\mathrm{5}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$, ${\mathit{A}}_{\mathrm{217}}^{\mathrm{PS}}\mathrm{=}\mathrm{60}\mathrm{\pm }\mathrm{10}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$, ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}\mathrm{=}\mathrm{32}\mathrm{\pm }\mathrm{8}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$, ${\mathit{A}}_{\mathrm{143}}^{\mathrm{PS}}\mathrm{=}\mathrm{75}\mathrm{\pm }\mathrm{8}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$, and ${\mathit{A}}_{\mathrm{100}}^{\mathrm{PS}}\mathrm{=}\mathrm{220}\mathrm{\pm }\mathrm{53}\hspace{0.17em}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}$. These estimates of the Poisson power are in good agreement with the predictions given in Sect. 7.3.1 for the 100 and 143 GHz channels. In the latter, radio sources below the Planck flux cuts dominate the Poisson power, which can be reliably estimated from existing source counts measurements.

We also consider modifying our model for the clustered part of the CIB. There have been a wealth of CIB models (e.g., Knox et al. 2001; Amblard & Cooray 2007; Hall et al. 2010; Pénin et al. 2012; Addison et al. 2013; Xia et al. 2012), most assuming that the dust is a biased tracer of the dark matter distribution, but differing in their parametrisation of the dust emissivity and its evolution, and their treatment of the dark matter power spectrum. Recent papers (Planck Collaboration XVIII 2011; Addison et al. 2013) have shown that the addition of Planck CIB measurements, when combined with other small scale probes including SPT, ACT, BLAST, and Herschel, rule out models that assume the underlying dark matter power spectrum is linear.

We therefore test a set of models that try to simultaneously fit the non-linear spectrum with one or more template spectra. We consider fixing the scale dependence to a power law, either ℓ^0.8 or ℓ^0.6, or extending the power law model to have a running of the index. The ℓ^0.8 model has been used in Shirokoff et al. (2011); Reichardt et al. (2012); Addison et al. (2012a); Dunkley et al. (2013), while the ℓ^0.6 more closely matches the CIB model of Addison et al. (2013). We find that a simple power law does not fit both the Planck and high-ℓ data sufficiently well, but that allowing the additional freedom of a running spectrum opens up the parameter space too much, with little improvement in goodness of fit, motivating our use of the varying γ^CIB model.

We also test a simple model using just a linear theory dark matter power spectrum, assuming that the non-linear power can be absorbed into the Poisson term. This results in an estimate of the Poisson level at 217 GHz that is inconsistent with ACT and SPT, so we do not use this model. For all these models we test the effect on cosmological parameters, using the Hubble constant as a test case, and find the effect on parameters to be small. This is also investigated in Planck Collaboration XVI (2014).

7.4. Consistency of the Planck 70 GHz data

Figure 11 shows individual frequency spectra from 70 to 353 GHz. We only use data from the 100 to 217 GHz channels to form the high-ℓ likelihoods, but here we compare cosmological parameters derived from the Planck 70 GHz channel alone. Maps at 70 GHz are easier to characterise than the higher frequency channels in terms of instrumental and foreground properties, but the resolution and sensitivity are lower. The 70 GHz noise properties are in general well described by a simple three parameter model involving 1 /f and white noise contributions (Planck Collaboration II 2014). The 70 GHz channel also has the least diffuse foreground emission at large scales (Planck Collaboration XII 2014), and the extragalactic source contribution is dominated by radio galaxies whose emission is well known at these frequencies (Planck Collaboration XXVIII 2014). We adopt a Galactic plane cut leaving ~70% of the sky for the analysis (CS70 from Planck Collaboration XII 2014) to which we add a point source mask optimised for 70 GHz. In Appendix D.2 we describe how cosmological parameters are estimated from the 70 GHz channel, which are summarized in Fig. 24. Accounting for the lower sensitivity and angular resolution at 70 GHz, which translate into a narrower multipole range (ℓ < 1200), the parameter distributions are consistent with the reference values.

7.5. Consistency with power spectra of CMB maps obtained by component separation methods

The likelihoods we consider in this paper account for component separation by modelling the multi-frequency data at the power spectrum level, to fully exploit the signal at the smallest scales probed by Planck. We can compare the results to those derived from an alternative approach, measuring the power spectrum of CMB maps estimated from component separation techniques. Here we present results obtained with four CMB maps, derived using methods referred to as Commander-Ruler, SMICA, NILC, and SEVEM, described in detailed in the accompanying paper (Planck Collaboration XII 2014). These maps are weighted with the union of the confidence masks of the different component separation methods, covering 72% of the sky (mask U73 of Planck Collaboration XII 2014), with a further cosine apodisation of 30′. We compare their angular power spectra and cosmological parameters with those from CamSpec. To estimate the power spectra we use the XFaster method, an approximation to the iterative, maximum likelihood, quadratic band power estimator based on a diagonal approximation to the quadratic Fisher matrix estimator (Rocha et al. 2011, 2010). The noise bias is estimated using difference maps, as described in Planck Collaboration XII (2014). The resulting spectra are shown in Fig. 29, and agree well out to scales ℓ_max = 2000, even though the agreement is less striking for the Commander-Ruler small-scale spectrum.

We then estimate cosmological parameters using a Gaussian correlated likelihood derived from these band-powers. To model the residual extragalactic foregrounds, we consider two nuisance parameters: A_ps, the amplitude of a Poisson component, and A_Cℓ, the amplitude of a clustered component, scaling a term with shape D_ℓ ∝ ℓ^0.8. Figure 30 compares the parameters obtained as a function of ℓ_max for each method, compared to the values from the high-ℓ likelihoods. The results are consistent to 1σ or better for ℓ_max = 1500 and ℓ_max = 2000. Despite adopting a simple two-parameter model for the extragalactic foregrounds, the likelihood using a CMB map as input data works reasonably well, and may be further exploited with analysis of simulations, improved extragalactic foreground modelling, and the development of an error model.

7.6. Consistency with high-ℓ polarisation

Our final consistency test is illustrated in Fig. 31, showing the polarisation power spectra derived from the Planck data. Both the TE and EE cross-spectra are shown, in bins of width of Δℓ = 40. These spectra are computed by performing a uniformly weighted average of all detector sets combinations at 70 × 100, 100 × 143, 100 × 217, and 143 × 217 GHz. We use the temperature beam window functions for beam deconvolution. For the analysis, we applied CL39 to the temperature maps, and discarded 60% of the polarisation data (i.e., f_sky = 0.4). Other than masking, no efforts have been made to subtract foreground contributions or take instrumental effects such as leakage into account. Despite the substantial masking applied, we see evidence for residual contributions of non-cosmological origin. Besides demonstrating the potential of Planck to deliver high quality polarisation maps and spectra (with the limitations explained in Planck Collaboration VI 2014), the figure demonstrates the high level of consistency of these polarisation spectra between themselves, and with the prediction from the model fit using just the temperature spectrum shown in Fig. 15.

As discussed in Planck Collaboration VI (2014) and Planck Collaboration II (2014), at present, the HFI and LFI polarisation spectra at low multipoles are affected by systematic errors that cause biases, which will need to be accurately modelled or removed for the next Planck release. However, these systematics rapidly become unimportant at higher multipoles. While not yet fit for cosmological parameter analysis, the consistency at the level of a few μK of these Planck polarisation spectra adds to our confidence in the analysis of temperature data. It shows that within the ΛCDM framework, the cosmological parameters estimated from Planck temperature data are not strongly affected by the uncertainties in the modelling of unresolved foregrounds.

8. Low-ℓ likelihood

At low multipoles (ℓ ≲ 50), the distribution of the estimated C_ℓs is not well approximated by a Gaussian due to the limited degrees of freedom per ℓ (e.g., Efstathiou 2004). However, both the CMB signal, s, and instrumental noise, n, are individually nearly Gaussian distributed at the map level, provided that foreground emission and instrumental systematics effects are negligible (e.g., Planck Collaboration XXIII 2014; Planck Collaboration XXIV 2014), and the actually observed map, m = s + n, is therefore also nearly Gaussian distributed. Under this assumption, the CMB power spectrum likelihood is given by $\mathrm{ℒ}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right)}\mathrm{=}\mathit{P}\mathrm{\left(}{m}\mathrm{\right|}{\mathit{C}}_{\mathit{\ell }}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{\pi }}^{\mathit{n}\mathit{/}\mathrm{2}}\mathrm{|}{M}{\mathrm{|}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{{m}}^{\mathrm{t}}\hspace{0.17em}{{M}}^{-1}{m}\right)\mathit{,}$(21)where n is the number of observed pixels, M(C_ℓ) = C(C_ℓ) + N is the data covariance matrix, and C and N are the CMB and noise covariance matrices, respectively.

In the general case, the data vector m includes both temperature (T) and linear polarisation (Q, U) Stokes parameter maps. Pixels exhibiting high foreground contamination are removed by masking, such that the data vector is restricted to the subset of valid pixels, m = (T_i1,T_i2,...,T_nT,Q_j1,Q_j2,...Q_nP,U_j1,U_j2,...U_nP). The corresponding rows and columns are removed from M, effectively corresponding to marginalising over the masked region of the sky. Note that in general, n_T ≠ n_P, and the sets of indexes of temperature and polarisation measurements are different. We assume the same number of pixels in Q and U, although this is not a requirement.

The signal covariance matrix can be written symbolically as ${C}\mathrm{=}\left(\begin{array}{c}\\ {⟨\mathit{TT}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{T}}\mathrm{\times }{\mathit{n}}_{\mathit{T}}\mathrm{)}}& {⟨\mathit{TQ}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{T}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}& {⟨\mathit{TU}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}\\ {⟨\mathit{QT}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{T}}\mathrm{)}}& {⟨\mathit{QQ}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}& {⟨\mathit{QU}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}\\ {⟨\mathit{UT}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{T}}\mathrm{)}}& {⟨\mathit{UQ}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}& {⟨\mathit{UU}⟩}_{\mathrm{(}{\mathit{n}}_{\mathit{P}}\mathrm{\times }{\mathit{n}}_{\mathit{P}}\mathrm{)}}\end{array}\right)\mathit{,}$(22)where the signal correlations for the temperature component are explicitly given by $\mathrm{⟨}{\mathit{T}}_{{\mathit{i}}_{\mathrm{1}}}{\mathit{T}}_{{\mathit{i}}_{\mathrm{2}}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{\ell }\mathrm{=}\mathrm{2}}^{{\mathit{\ell }}_{\max }}\frac{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathit{Ĉ}}_{\mathit{\ell }}{\mathit{P}}_{\mathit{\ell }}\mathrm{\left(}{\mathit{\theta }}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}}\mathrm{\right)}\mathrm{+}{{N}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}}\mathit{.}$(23)Here P_ℓ are the Legendre polynomials, and θ_i1i2 is the angle between the centres of pixels i₁ and i₂. Similar expressions are available for the polarisation correlations (e.g., Tegmark & de Oliveira-Costa 2001). The effect of the (azimuthally symmetric) instrumental beam, b_ℓ, and pixel window function, w_ℓ, are encoded in ${\mathit{Ĉ}}_{\mathit{\ell }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{th}}{\mathit{b}}_{\mathit{\ell }}^{\mathrm{2}}{\mathit{w}}_{\mathit{\ell }}^{\mathrm{2}}$.

The main problem with the likelihood expression given in Eq. (21) is its high computational cost. This is determined by the matrix inversion and determinant evaluations, both of which scale as with N = n_T + 2n_P. In practice, this approach is therefore limited to coarse pixelisations, N_side ≤ 16, which reliably only supports multipoles below ℓ ≲ 30. On the other hand, the Gaussian approximation adopted by the high-ℓ likelihood is not sufficiently accurate for the stringent requirements of Planck below ℓ ≲ 50. In the next section, we therefore describe a faster low-ℓ likelihood estimator, based on Gibbs/MCMC sampling, which allows us to exploit the full range up to ℓ ≤ 50 with low computational cost, while additionally supporting physically motivated foreground marginalisation.

Page et al. (2007) pointed out that the temperature and polarisation parts of the likelihood can be separated and evaluated independently, under the assumption of negligible noise in temperature and in the temperature-polarisation cross correlations (i.e., the TQ and TU blocks of the pixel level noise covariance matrices). Further assuming vanishing primordial B modes and TB correlations, the TE correlations can be accounted for by redefining the modified Q and U maps as $\begin{array}{ccc}\mathit{Q}& \mathrm{\to }& \\ \mathit{U}& \mathrm{\to }& \end{array}$where ${}_{\mathrm{\pm }\mathrm{2}}\mathit{Y}_{\mathit{\ell m}}$ are spin weighted spherical harmonics and ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}$ are the harmonic coefficients of the signal in the temperature map. One can show by direct substitution that these modified Q and U maps are free of temperature correlations. The polarisation likelihood can be then computed independently from the temperature likelihood and, possibly, at lower resolution to save computational expenses. We test this strategy in Sect. 8.2, and adopt it for the current release of the Planck likelihood.

8.1. Low-ℓ temperature likelihood

As discussed above, we do not implement the likelihood expression given in Eq. (21) directly, due to its high computational cost and limited flexibility with respect to foreground modelling. Instead, we adopt the Gibbs sampling approach (Eriksen et al. 2004; Jewell et al. 2004; Wandelt et al. 2004), as implemented by the Commander code (Eriksen et al. 2008), which allows both for physically motivated component separation and accurate likelihood estimation. A similar Gibbs sampling method was used to estimate the low-ℓ temperature likelihood for WMAP (Dunkley et al. 2009; Larson et al. 2011), although not simultaneously accounting for component separation.

8.1.1. Methodology

We start by generalizing the above data model to include both multi-frequency observations and a set of foreground signal terms, ${{d}}_{\mathit{\nu }}\mathrm{=}{s}\mathrm{+}\sum _{\mathit{i}}{{f}}_{\mathit{\nu }}^{\mathit{i}}\mathrm{+}{{n}}_{\mathit{\nu }}\mathit{.}$(26)Here d_ν denotes the observed sky map at frequency ν, and ${{f}}_{\mathit{\nu }}^{\mathit{i}}$ denotes a specific foreground signal component. As above, the CMB field is assumed to be a Gaussian random field with power spectrum C_ℓ, and the noise is assumed Gaussian with covariance matrix N_ν. The foreground model can be adjusted as needed for a given data set, and a full description of the model relevant for Planck is presented in Planck Collaboration XII (2014). In short, this consists of a single low-frequency foreground component (i.e., the sum of synchrotron, anomalous microwave emission, and free-free emission), a carbon monoxide (CO) component, and a thermal dust component, in addition to unknown monopole and dipole components at each frequency.

Given this data model, we map out the full posterior distribution, P(s,fⁱ,C_ℓ | d), using a Monte Carlo sampling algorithm called Gibbs sampling. Directly drawing samples from P(s,fⁱ,C_ℓ | d) is computationally prohibitive, but this algorithm achieves the same by iteratively sampling from each corresponding conditional distribution,

$\begin{array}{ccc}& & {s}\mathrm{\leftarrow }\mathit{P}\mathrm{\left(}{s}\mathrm{\right|}{f}\mathit{,}{\mathit{C}}_{\mathit{\ell }}\mathit{,}{d}\mathrm{)}\\ & & {f}\mathrm{\leftarrow }\mathit{P}\mathrm{\left(}{f}\mathrm{\right|}{s}\mathit{,}{\mathit{C}}_{\mathit{\ell }}\mathit{,}{d}\mathrm{)}\\ & & {\mathit{C}}_{\mathit{\ell }}\mathrm{\leftarrow }\mathit{P}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right|}{s}\mathit{,}{{f}}^{\mathit{i}}\mathit{,}{d}\mathrm{)}\mathit{.}\end{array}$It is straightforward to show that P(s | f,C_ℓ,d) is a multivariate Gaussian distribution, and P(C_ℓ | s,fⁱ,d) is an inverse Gamma distribution. The foreground distribution, P(f | s,C_ℓ,d), does not have a closed analytic form, but can easily be mapped out numerically (Eriksen et al. 2008). Thus, all three distributions are associated with simple textbook sampling algorithms.

For CMB power spectrum and likelihood estimation using Gibbs sampling, the crucial intermediate product from the above sampling process is the ensemble of CMB sky samples, s^k. Each individual sample corresponds to one possible CMB realisation consistent with the observed data. In the absence of sky cuts, foreground contamination and instrumental noise, this map is identical to the true sky. In that case, the likelihood as a function of C_ℓ is determined by cosmic variance alone, and given by an inverse gamma distribution, ${\mathrm{ℒ}}^{\mathit{k}}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right)}\mathrm{\propto }\frac{{\mathit{\sigma }}_{\mathit{\ell ,k}}^{\frac{\mathrm{2}\mathit{\ell }\hspace{0.17em}\mathrm{-}\hspace{0.17em}\mathrm{1}}{\mathrm{2}}}}{{\mathit{C}}_{\mathit{\ell }}^{\frac{\mathrm{2}\mathit{\ell }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}{\mathrm{2}}}}{\mathrm{e}}^{\mathrm{-}\frac{\mathrm{2}\mathit{\ell }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}{\mathrm{2}}\frac{{\mathit{\sigma }}_{\mathit{\ell }}\mathit{,k}}{{\mathit{C}}_{\mathit{\ell }}}}\mathit{.}$(27)Here we have introduced the realisation specific power spectrum, ${\mathit{\sigma }}_{\mathit{\ell ,k}}\mathrm{\equiv }\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}{\sum }_{\mathit{\ell }\mathrm{=}\mathrm{-}\mathit{m}}^{\mathit{m}}\mathrm{|}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{k}}{\mathrm{|}}^{\mathrm{2}}$, where ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{k}}$ are the spherical harmonic coefficients of s^k. In the case of realistic data, we need to marginalise over uncertainties due to sky cuts, foregrounds, and instrumental noise. Hence, Eq. (27) is replaced by an average over all possible sampled CMB realisations, $\begin{array}{ccc}\mathrm{ℒ}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right)}\mathrm{\propto }\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{samp}}}{\mathrm{ℒ}}^{\mathit{k}}\mathrm{\left(}{\mathit{C}}_{\mathit{\ell }}\mathrm{\right)}\mathit{.}& & \end{array}$(28)This expression is known as the Blackwell-Rao estimator (Chu et al. 2005), and is guaranteed to converge to the exact likelihood as the number of samples, N_samp, increases. Note that the normalization factor in this expression is unknown, but since the likelihood function is only used to compare different models through an effective likelihood ratio, this factor is irrelevant for actual calculations.

8.1.2. Data selection and preprocessing

As described in Planck Collaboration XII (2014), we include Planck frequencies between 30 and 353 GHz in the low-ℓ likelihood. Each frequency map is downgraded from its native resolution to a common resolution of 40′, and projected onto an N_side = 256 HEALPix grid. Uncorrelated Gaussian regularization noise is added to each frequency channel map, with an RMS proportional to the spatial mean of the instrumental noise of the corresponding channel, ⟨σ_ν⟩, conserving relative signal-to-noise between channels. The regularization noise level at frequency ν is $\mathrm{5}\mathit{\mu }\mathrm{K}\mathrm{·}⟨{\mathit{\sigma }}_{\mathit{\nu }}⟩\mathit{/}{}^{\mathrm{⟨}}\mathit{\sigma }_{\mathrm{143}\mathrm{GHz}}{}^{\mathrm{⟩}}$. The purpose of this is to make the results insensitive to unmodelled features at scales comparable to and beyond the smoothing scale of 40′, in addition to improve the convergence speed of the Gibbs sampler. The resulting signal-to-noise is unity at ℓ ~ 400, and the additional uncertainty due to the regularization noise is less than 0.2 μK² below ℓ = 50, and less than 1 μK² below ℓ = 100.

To study the stability of the low-ℓ likelihood with respect to sky fraction, we constructed a suite of five different masks, leaving between 81 and 100% of the sky for the analysis, and for completeness we include the WMAP KQ85 mask (75% of the sky), as a sixth case. These low-ℓ masks are distinct from those employed for the high-ℓ likelihood, and are produced in a preliminary Commander full-sky analysis in which we estimate individual foreground components and residual χ² values per pixel. These maps are thresholded at various levels to produce a useful range of sky fractions.

For each mask, we fit a two-parameter amplitude and tilt power spectrum model of the form ${\mathit{C}}_{\mathit{\ell }}\mathrm{\left(}\mathit{q,n}\mathrm{\right)}\mathrm{=}\mathit{q}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{fid}}{{}^{\mathrm{(}}\mathit{\ell }\mathit{/}{\mathit{\ell }}_{\mathrm{0}}{}^{\mathrm{)}}}^{\mathit{n}}$, using the low-ℓ likelihood between ℓ_min = 2 and ℓ_max = 49, where ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{fid}}$ is the best-fit PlanckΛCDM spectrum, and ℓ₀ = (ℓ_min + ℓ_max) / 2. The resulting distributions are shown in Fig. 32 for three masks, covering 77.4 (A; WMAP KQ85), 87.5 (B) and 97.2% (C) of the sky, respectively. The internal agreement is excellent, with parameters differing by less than 0.3σ between the very aggressive mask A and the conservative mask C. While any of these masks would establish an acceptable likelihood, we adopt Mask B as our fiducial mask for two reasons. On the one hand, the parameter uncertainties obtained with Mask B are only 4% larger than those obtained for the minimal Mask A, indicating that both nearly saturate the cosmic variance limit. On the other hand, analysis of realistic simulations indicate the presence of statistically significant map residuals near the Galactic plane that are accepted by Mask A, but rejected by Mask B (Planck Collaboration XII 2014). The latter therefore represents a good compromise between rejecting foreground residuals and maximizing statistical power.

We include 100 000 Gibbs samples in the likelihood estimator, ensuring excellent convergence characteristics for the Blackwell-Rao estimator for ℓ < 50.

8.2. Low-ℓ polarisation likelihood

The present Planck data release includes only temperature data. In this release, we therefore supplement the Planck likelihood with the 9 year WMAP polarisation likelihood⁹ derived from the WMAP polarisation maps at 33, 41, and 61 GHz (Ka, Q, and V bands) (Page et al. 2007; Bennett et al. 2013). However, we introduce one modification to this pixel-based likelihood code, replacing the spherical harmonics coefficients of the temperature field, ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}$, in Eqs. (24) and (25) with those derived from the Planck temperature map derived by Commander, for which the Galactic plane has been replaced with a Gaussian constrained realisation.

In Fig. 33, we compare constraints on τ and A_s as derived with this split likelihood with those obtained through an exact brute-force evaluation of Eq. (21), simultaneously including temperature and polarisation measurements at N_side = 16. The two methods produce almost indistinguishable results.

In Appendix E we assess the robustness of the WMAP polarisation likelihood with respect to dust contamination, by replacing the WMAP polarised dust template with the far more sensitive HFI 353 GHz polarisation map. We find that the optical depth to reionisation, τ, is reduced by about 0.5 − 1σ, depending on the template removal method adopted. However, since the Planck polarisation maps are not included from the current data release, we adopt the WMAP polarisation likelihood without further changes for now, and will return to this topic in the next data release.

8.3. Low-ℓ power spectrum – consistency and robustness

In this section, we present the low-ℓPlanck CMB temperature power spectrum derived using the Commander approach described above, and assess its robustness through comparisons with three alternative foreground-cleaned Planck CMB maps (NILC, SEVEM, and SMICA; Planck Collaboration XII 2014), as well as with the 9 year WMAP ILC temperature map (Bennett et al. 2013).

As a first consistency test, we compute the power spectrum from each map using Bolpol (Gruppuso et al. 2009), an implementation of the quadratic maximum-likelihood power spectrum estimator (Tegmark 1997). Each map is smoothed to an effective resolution of 329.81′ FWHM, to suppress aliasing from high multipoles (Keskitalo et al. 2010), and repixelised on an N_side = 32 HEALPix grid. Gaussian white noise with a variance of 4 μK² is added to each map to regularise the noise covariance matrix.

Here we adopt the U78 common mask, defined in Planck Collaboration XII (2014), to exclude regions of high Galactic emission, leaving 78% of the sky for analysis. We remove the observed monopole and dipole in the mask. The resulting power spectra up to ℓ ≤ 64 are shown in the top panel of Fig. 34, while the bottom panel shows the power spectrum residuals of each map relative to the Commander map. Note that the same noise realisation was added to each map, and the regularization noise therefore contributes little in this plot. For the different internally-derived Planck maps, no residual spectrum exceeds ≲100μK² and is typically ≲50μK² at ℓ ≳ 10. The WMAP spectrum exhibits significantly larger residuals, and are typically of the order of ~100 μK² at ℓ ≳ 30.

Figure 35 shows the Planck and WMAP temperature power spectra derived directly from the respective likelihood code, while Fig. 36 shows the corresponding constraints on the two-parameter amplitude-tilt model employed in Sect. 8.1.2, including multipoles between ℓ = 2 and 30. Neglecting the minor differences in the masks adopted by the two codes, these power spectra and parameter constraints are largely dominated by cosmic variance, and one should therefore expect the two distributions to be almost identical. Instead, from Fig. 36 we see that the WMAP low-ℓ spectrum is 2.5–3% higher than the Planck spectrum. For a detailed discussion of this discrepancy, including a comparison at higher ℓ, see Planck Collaboration XXXI (2014). Here we only note that the effect is robust with respect to foreground removal and power spectrum evaluation algorithms, and also point out that the effect at low multipoles is too large to be explained by uncertainties in the Planck transfer functions (Planck Collaboration II 2014; Planck Collaboration VII 2014) or calibration (Planck Collaboration V 2014; Planck Collaboration VIII 2014). Also note that the amplitude of the low-ℓ spectrum relative to the Planck best-fit model, (q,n) = (1,0), derived including the full multipole range between 2 ≤ ℓ ≤ 2500, is somewhat low in Fig. 36, with a best-fit amplitude of q ~ 0.9. This observation is discussed and quantified in greater detail in Sect. 9.3.

9. The Planck CMB spectrum and likelihood

9.1. Hybridisation of low- and high-ℓ likelihoods

The high-ℓ and low-ℓ likelihoods introduced in Sects. 2 and 8 each describe only a part of the full Planck data set. To estimate cosmological parameters from all the angular scales probed by Planck, they must be combined into a single likelihood function that describes all multipoles from ℓ = 2 to 2500.

In principle, it is desirable to include as many multipoles as possible in the low-ℓ likelihood, since it captures the full non-Gaussian structure of the likelihood. The Gaussian approximation for the likelihood using pseudo-spectra also improves at higher multipole due to the increasing number of degrees of freedom (Efstathiou 2004). For Planck we adopt a transition multipole of ℓ_trans = 50, a compromise between obtaining robust convergence properties for the low-ℓ likelihood, and ensuring that the Gaussian approximation holds for the high-ℓ likelihood (Hamimeche & Lewis 2009).

To combine the likelihoods, we must account for the weak correlations between the low- and high-ℓ components. We consider three options:

• 1.

  Sharp transition: the low-ℓ likelihood ends at ℓ_max = 49; the high-ℓlikelihood starts at ℓ_min = 50; no correlations are accounted for.

• 2.

  Gap: the low-ℓ likelihood ends at ℓ_max = 32; the high-ℓ likelihood starts at ℓ_min = 50; no correlations are accounted for, but the gap is sufficiently wide that any correlations are negligible.

• 3.

  Overlap with correction: the low-ℓ likelihood ends at ℓ_max = 70; the high-ℓ likelihood starts at ℓ_max = 50; the double-counting of the overlap region is accounted for by subtracting from the log-likelihood a contribution only including 50 ≤ ℓ ≤ 70 as evaluated by the Commander estimator. Under the assumption that no correlations extend from ℓ ≤ 50 to ℓ ≥ 70, this approach is exact).

We estimate cosmological parameters using all three methods, and find that the posterior means typically vary by <0.1σ. The largest variation is seen including the running of the spectral index of scalar perturbations, in which the posterior mean changes by 0.2σ. Further, all deviations at the 0.1–0.2σ level are seen for case 2 above, which excludes data compared to the other two; case 1 and 3 give nearly indistinguishable results. Since case 1 is implementationally simpler, and can be estimated more efficiently (see Sect. 8), we select this method, adopting a sharp transition at ℓ_max = 50.

9.2. The Planck power spectrum and ΛCDM constraints

Using the full Planck likelihood, we now present the final 2013 Planck CMB power spectrum. For this, we fix all nuisance parameters to their maximum-likelihood values. The resulting spectrum is shown in Fig. 37 together with the corresponding best-fit six-parameter ΛCDM model. The agreement between the observations and the model is excellent over most of the multipole range. Only at low ℓs is it possible to see a systematic offset in the form of a slight power deficit; this is addressed separately in the next section.

Table 8 provides a summary of the ΛCDM parameters derived using the methodology described in Planck Collaboration XVI (2014) from the Planck likelihood. Here we use the same prior ranges on all parameters as in Planck Collaboration XVI (2014). These are as in Table 5 and 6, except for ${\mathit{A}}_{\mathrm{143}}^{\mathrm{CIB}}$, A^tSZ, and A^kSZ, which are modified to [0,20], [0,10], and [0,10] respectively. Results are given for Planck alone, and in combination with the low-ℓ WMAP polarisation likelihood (Planck+WP). For each case, we report both posterior maximum and mean values. Uncertainties denote 68% confidence limits.

A detailed discussion of these results is presented in Planck Collaboration XVI (2014), including an analysis of extended cosmological models, and their compatibility with other astrophysical data sets. The bounds derived from Planck alone are significantly tighter than those from the 9 year WMAP data alone, and comparable or better than those inferred from WMAP combined with SPT and ACT observations. These new constraints provide a precision test of the ΛCDM model. In general, we find good agreement with results derived from other astrophysical data sets, although there are a few exceptions that are in moderate tension with Planck (Planck Collaboration XVI 2014).

Considering each of the six ΛCDM parameters in turn, we first note that Planck constrains the physical baryon density to Ω_bh² = 0.02207 ± 0.00033, which is in remarkable agreement with standard BBN predictions based on a determination of the primordial abundance of deuterium, Ω_bh² = 0.021 ± 0.001 (Iocco et al. 2009), but with a fractional uncertainty of 1.5 %, three times smaller than the BBN uncertainty. The physical density of dark matter is measured with a fractional uncertainty of 2.6%, providing new constraints on specific dark matter production scenarios. The single most precise parameter, however, is the angular size of the sound horizon at the last-scattering surface, θ_MC, which is measured with a fractional uncertainty of 0.065 % by Planck, improving on the combined WMAP, ACT, SPT, and SNLS constraint by a factor of two.

Next, given that no polarisation data are included in the current data release, it is remarkable that Planck alone constrains the optical depth to reionisation, τ, with a fractional error of 40%. This is made possible by Planck’s high angular resolution and sensitivity, which allows a high signal-to-noise measurement of lensing in the small-scale CMB power spectrum. This in turn breaks the well-known e${}^{\mathrm{-}\mathrm{2}\mathit{\tau }}\mathit{A}_{\mathrm{s}}$ degeneracy between τ and the amplitude of scalar perturbations, A_s. The fractional uncertainty on A_s from Planck alone is 7%.

Having sufficient power to measure τ from small angular scale temperature data, Planck naturally also provides very strong constraints on the spectral index of scalar perturbations, n_s, leading to a fractional uncertainty of 0.97 %. The scale-invariant Harrison-Zeldovich spectrum, n_s = 1, is ruled out at a significance of 4.1σ from the Planck temperature spectrum alone. The significance of this exclusion rises to 5.4σ when the low-ℓ WMAP polarisation is included. The analyses presented in Planck Collaboration XVI (2014) and Planck Collaboration XXII (2014) show that the preference for a (red) tilted primordial spectrum remains very strong also within most extensions beyond the minimal ΛCDM model. The implications of this results for inflationary models are discussed in Planck Collaboration XXII (2014).

With our choice of cosmological parameters, the Hubble parameter, H₀, and the fractional density of the cosmological constant, Ω_Λ = 1 − Ω_m, are derived parameters. They are probed by CMB observations mainly through their impact on θ_MC, and, to lesser extent, by the impact of Ω_Λ on the late-time integrated Sachs-Wolfe effect. Since θ_MC is accurately measured, a particular combination of H₀ and Ω_Λ is very well constrained by Planck, although in a model-dependent way; θ_MC depend on other cosmological parameters, such as the spatial curvature radius, neutrino masses, the number of relativistic degrees of freedom, or a possible dark energy equation of state parameter.

The results reported in Table 8 rely on the assumption of a flat ΛCDM cosmology with three neutrino species, two of which are assumed to massless and one featuring a small mass m_ν = 0.06 eV, reflecting the lower bound on neutrino masses imposed by neutrino oscillation experiments. Under these assumptions, Planck finds preferred ranges for H₀ and Ω_Λ that are lower than previous CMB experiments. For instance, Planck+WP gives H₀ = 67.3 ± 1.2 km s^-1Mpc^-1, to be compared with 70.5 ± 1.6 km s^-1Mpc^-1 for the combined WMAP9+eCMB data set presented by Hinshaw et al. (2013). The underlying cosmology in the two analyses is the same, excepted for the small neutrino mass introduced in our default ΛCDM model. However, if we assume all three neutrino species to be massless, our best-fit and mean values for H₀ increase only by 0.6 km s^-1Mpc^-1. Thus, the tension is clearly driven by the data rather than by theoretical assumptions. Planck Collaboration XVI (2014) shows that our results for H₀ and Ω_Λ are in very good agreement with Baryon Acoustic Oscillation data, but in moderate tension with other cosmological probes. For instance, our Planck+WP bounds on H₀ disagree at the 2.5σ level with direct determinations of the Hubble parameter using Cepheid and supernovae (Riess et al. 1998) or quasar time delays (Suyu et al. 2013), as well as with the results of the Carnegie Hubble Program (Freedman et al. 2012). Our bounds on Ω_Λ are in a slight 2σ tension with the results of the SNLS supernovae collaboration (Conley et al. 2011; Sullivan et al. 2011), although in better agreement with the Union2.1 compilation (Suzuki et al. 2012). Our combined determination of σ₈ and Ω_m shows larger tension with recent data based on cosmic shear or cluster count techniques. On the other hand the Planck best-fit ΛCDM model is in good agreement with the halo power spectrum derived from the luminous red galaxy catalogue of the Sloan Digital Sky Survey (Reid et al. 2010), especially when the analysis is restricted to linear scales.

9.3. Significance of the low-ℓ tension with ΛCDM models

From the above discussion, it is clear that the ΛCDM framework provides an excellent model for most of the Planck data. However, as noted in Sect. 8.3 and seen in Fig. 37, the low-ℓPlanck temperature power spectrum appears to be in some tension with the best-fit PlanckΛCDM model, which for Planck is almost exclusively determined by the small-scale spectrum. In this section we assess the significance and impact of this tension between low and high ℓs using three different statistical tests.

We start by applying a modified Hausman test (Polenta et al. 2005; Planck Collaboration II 2014) to the low-ℓ spectra derived from the four foreground-cleaned Planck maps (Planck Collaboration XII 2014) and the 9-year WMAP ILC map, using multipoles up to ℓ_max = 32. This test uses the statistic s₁ = sup_rB(ℓ_max,r), where

$\begin{array}{ccc}& & \\ & & {\mathit{H}}_{\mathit{\ell }}\mathrm{=}\frac{\widehat{{\mathit{C}}_{\mathit{\ell }}}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}}{\sqrt{\mathrm{Var}\hspace{0.17em}\widehat{{\mathit{C}}_{\mathit{\ell }}}}}\mathit{,}\end{array}$and Ĉ_ℓ and C_ℓ denote the observed and model power spectra, respectively. Intuitively, this statistic measures the relative bias between the observed spectrum and model, measured in units of standard deviations, while taking the so-called “look-elsewhere effect” by maximizing s₁ over multipole ranges into account. We use realistic Planck “FFP6” simulations (Planck Collaboration I 2014) to derive the empirical distribution of s₁ under the null hypothesis. Figure 38 compares the results obtained from the data with the simulation distribution, and Table 9 lists significances. As measured by this statistic, we see that a negative bias is found in the low-ℓPlanck power spectrum relative to the ΛCDM model at the 99% confidence level.

For the WMAP ILC map the significance of the negative bias nominally decreases to 93%. This is consistent with the findings in Sect. 8.3, where it was shown that the WMAP temperature power spectrum is 2.5–3 % higher than the Planck spectrum at low ℓ’s. However, as discussed in Planck Collaboration XXXI (2014), a similar amplitude difference between the two experiments is also seen at smaller scales. Since the current test compares the observed WMAP data with the best-fit PlanckΛCDM model, the present test is not optimal for assessing internal consistency between low and high ℓs within the WMAP data.

Next, to obtain a quantitative measure of the relative power discrepancy between low and high ℓs, we fit the two-parameter amplitude–tilt power spectrum model (see Sect. 8.1.2) to the Planck data using the low-ℓ likelihood restricted to various multipole ranges defined by 2 ≤ ℓ ≤ ℓ_max, where ℓ_max is allowed to vary. Thus, this measures the amplitude of the low-ℓ spectrum relative to the best-fit PlanckΛCDM spectrum, which is driven by the smaller angular scales. Figure 39 shows the resulting constraints on the power spectrum amplitude, q, as a function of ℓ_max, after marginalising over the tilt, n. For comparison, we also show similar constraints derived using the low-ℓ WMAP temperature likelihood up to ℓ = 30. The best-fit amplitude is q ~ 0.9 for ℓ_max = 20–35, different from unity at a statistical significance of 2–2.5σ by this measure. The WMAP spectrum shows a consistent behaviour, up to the same overall scaling factor of 2.5–3% between Planck and WMAP discussed above. We have verified that these results are insensitive to the (well-known) low quadrupole moment by excluding C₂ from the analysis; the large cosmic variance of this particular mode results in a low overall statistical weight in the fit.

Finally, we assess the impact of the low-ℓ power deficit on the ΛCDM model estimated using the Planck likelihood¹⁰ (augmented with the WMAP polarisation likelihood). We fit a low-ℓ rescaling amplitude, A_low for ℓ <ℓ_low jointly with the ΛCDM parameters, i.e., ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}{\mathit{A}}_{\mathrm{low}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{\Lambda }\mathrm{CDM}}$ for ℓ <ℓ_low and ${\mathit{C}}_{\mathit{\ell }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{\Lambda }\mathrm{CDM}}$ for ℓ ≥ ℓ_low. Figure 40 shows the resulting posterior distributions for A_low for ℓ_low = 32 (green) and ℓ_low = 49 (blue). The purple line shows the same when replacing the Planck low-ℓ likelihood with the WMAP low-ℓ likelihood (ℓ_low = 32). The corresponding best-fit values are A_low = 0.899 ± 0.046 (Planck; ℓ_low = 32), A_low = 0.953 ± 0.033 (Planck; ℓ_low = 49) and A_low = 0.953 ± 0.048 (WMAP; ℓ_low = 32), respectively. As already noted in Sect. 8.3, the observed deviations from unity are too large to be explained by the <1% uncertainties in the Planck transfer functions (Planck Collaboration II 2014; Planck Collaboration VI 2014).

In Fig. 41 we show the posterior distributions for Ω_ch², n_s and H₀ after marginalising over A_low for ℓ_low = 49. (Adopting ℓ_low = 32 results in negligible differences for all parameters except A_low). Shifts of 0.6–1σ are observed compared to the reference model, A_low = 1. We note that H₀, which already has a “low” value (for a detailed discussion, see Planck Collaboration XVI 2014), prefers an even lower value when allowing a rescaling of the low-ℓ spectrum. As a final test, we replace the entire low-ℓ likelihood, both temperature and polarisation, with a Gaussian prior on the optical depth of reionisation, τ = 0.089 ± 0.014, matching the WMAP measurement (Hinshaw et al. 2013). The resulting posteriors are shown as purple lines in Fig. 41, and agree well with the case including a low-ℓ scaling factor, but are, in fact, slightly further away from the reference model. Although not very significant in an absolute sense, these results do indicate that the high-ℓ likelihood is challenged in finding models that also fit the low-ℓ power spectrum.

To summarize, we have phenomenologically quantified a tension between the low-ℓ CMB power spectrum at ℓ ≲ 40 and the Planck best-fit ΛCDM model. Its significance varies between 2.5 and 3σ depending on the estimator used. The effect is seen in all four Planck foreground-cleaned CMB maps with little variation. It is also present in the 9-year WMAP data, although an overall amplitude difference of 2.5–3% between the data sets complicate a direct comparison. To make further progress, one would seek to establish a physical model that predicts a low-ℓ power deficit compared to high ℓ’s, and that may also predict other observable effects, which may be tested with cosmological data. Such a model may be related to the tentative detections of violations of statistical isotropy discussed in Planck Collaboration XXIII (2014), e.g., the low CMB temperature variance, the hemispherical power asymmetry, or the alignment between the quadrupole and octopole moments.

10. Discussion and conclusions

We have presented the Planck likelihood, which provides a detailed and accurate characterisation of the two-point statistics of the CMB temperature field, accounting for all significant sources of uncertainty; statistical, instrumental, and astrophysical. This likelihood function allows us to present an estimate of the CMB temperature power spectrum that spans more than three decades in ℓ with unprecedented precision; a spectrum that saturates the cosmic variance limit at all scales ≳0.1°, nearly exhausting the information content of the temperature anisotropies, and, in fact, becoming limited by uncertainties due to astrophysical foreground modelling. This is precisely what was originally promised at the time when Planck was selected by ESA in March 1996.

On large angular scales, ℓ < 50, the Planck likelihood is based on a Gibbs sampling approach that allows joint CMB power spectrum and component separation analysis, while accurately marginalising over a physically motivated foreground model constrained by the 30–353 GHz Planck frequencies. On intermediate and small scales, the Planck likelihood employs a fine-grained set of cross-spectrum combinations among the 100, 143, and 217 GHz detector maps to constrain the high-ℓ CMB power spectrum, ensuring that no noise bias can compromise the results, while at the same time allowing for physical foreground modelling in terms of power spectrum templates. This emphasis on physical foreground modelling has made it possible to combine the full power of the Planck data with observations from higher-ℓ CMB experiments.

We have validated our results through an extensive suite of consistency and robustness analyses, propagating both instrumental and astrophysical uncertainties to final parameter estimates. Further, we have studied in detail the well-known degeneracies that exist between the foreground and cosmological parameters at high ℓs when only including Planck observations, and shown that they have only a weak impact on cosmological conclusions.

On a more detailed level, we draw the following conclusions:

• The consistency between power spectra measured within thedifferent frequency channels is remarkable. In thesignal-dominated regime for single detectors, at ℓ ≲ 1000, thecross-spectra show an RMS dispersion of a few μK² in multipolebands of δℓ = 31. This confirms the relative calibration of the 100, 143,and 217 GHz detectors to ~0.2%.

• The differences, , between the 143 × 143, 143 × 217, and 217 × 217 cross-spectra averaged over multipole bands of δℓ ≈ 31 have a dispersion over 800 ≤ ℓ ≤ 1500 of 9, 5, and 5 μK², respectively, after subtracting the best-fit foreground model. This dispersion is not primarily of instrumental origin, but can be predicted from a model of the chance correlations between foregrounds and CMB fluctuations.

• At high ℓs, the power spectrum of the four foreground-cleaned CMB maps derived through component separation are consistent within their uncertainties. The cosmological parameters derived from these maps are consistent with those estimated by the Planck likelihood for ℓ ≲ 2000, despite very different foreground models.

• At low ℓs, the power spectrum differences among the four foreground-cleaned CMB maps are below 50 μK² for nearly every single multipole. Residuals with respect to the 9 year WMAP ILC map are slightly larger, typically 100 μK² or more. A detailed comparison between Planck and WMAP reveals a systematic power spectrum amplitude difference at the 2–3% level that cannot be accounted for within the Planck instrumental error budget. This is consistent with the findings presented in Planck Collaboration XXXI (2014).

• Parameters derived from the 70 GHz Planck frequency map are in excellent agreement with the reference results derived using the Planck likelihood; when the latter is limited to ℓ ≤ 1000, the agreement is even more striking. This confirms the strong internal consistency between the LFI and HFI instruments.

• The best-fit ΛCDM model derived from the Planck likelihood predicts TE and EE spectra in good agreement with the measured polarisation signature over a broad range of frequencies (70 to 217 GHz) and multipoles (ℓ ≲ 1000). At 100, 143, and 217 GHz, the instrumental noise in the EE spectrum is at the μK² level for ℓ ≲ 1000, and the visible (small) differences between the spectra suggests different levels of foreground contribution.

• We report a tension between the Planck best-fit ΛCDM model and the low-ℓ spectrum in the form of a power deficit of 5–10% at ℓ ≲ 40, with a statistical significance of 2.5–3σ. Thus, while the minimal ΛCDM model provides an outstanding fit for intermediate and small angular scales, this tension may suggest that the model is incomplete. In this respect, it is worth noting that other, but possibly related, anomalies have been reported in a companion paper studying statistical isotropy in the Planck sky maps at statistically significant levels.

In summary, we find that the majority of the Planck data can be described by a minimal six-parameter ΛCDM model with a very high degree of accuracy. Within this model the statistical uncertainties are dominated by astrophysical foreground modelling by scales of ℓ ≃ 1500. At lower ℓs, the unprecedented quality of the Planck data is such that the only fundamental limit is that we can only observe one CMB sky. In other words, Planck is cosmic variance dominated at ℓ ≲ 1500, extragalactic foreground dominated at ℓ ≳ 1500, and dominated nowhere by instrumental noise or systematic errors.

Using only Planck temperature data, we report a detection of n_s< 1 at more than 4σ confidence, which by itself is already a bit stronger than the limit derived from WMAP, SPT, ACT, and SNLS combined. Complementing the Planck observations with the 9 year WMAP polarisation data increases the significance further to 5.4σ. The multipole range above ℓ > 1500 is crucial for constraining possible extensions to the minimal ΛCDM model; for a detailed exploration of a wide range of such models, see Planck Collaboration XVI (2014) and Planck Collaboration XXII (2014). There we report some tensions among the CMB damping tail parameters, including Ω_K, n_run, and Y_P. However, none of these indicate significant departures from the ΛCDM framework.

In the near future, we will extend our analysis to produce a cosmic variance limited likelihood and power spectrum reaching to higher multipoles. To some extent, this will be achieved through more sophisticated astrophysical foreground modelling, and by exploiting additional frequency information. However, the two major steps forward will be, first, to include the Planck polarisation observations in the likelihood analysis, and, second, to exploit the full data set generated by the two Planck instruments. The amount of HFI data available for analysis is nearly double that, which is presented here, and the LFI instrument is still observing at the time of writing.

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (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.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. We acknowledge the use of the CLASS Boltzmann code (Lesgourgues 2011) and the Monte Python package (Audren et al. 2013) in earlier stages of this work. The likelihood code and some of the validation work was built on the library pmclib from the CosmoPMC package (Kilbinger et al. 2011). This research used resources of the IN2P3 Computer Center (http://cc.in2p3.fr) as well as of the Planck-HFI data processing centre infrastructures hosted at the Institut d’Astrophysique de Paris (France) and financially supported by CNES.

References

Appendix A: High-ℓ likelihood details

Appendix A.1: Power spectra and the coupling matrix

We denote the pixel weight function for temperature by ${\mathit{w}}_{\mathit{i}}^{\mathit{T}}$. The pseudo-spectra of Eq. (1) are constructed using the following: $\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}\mathrm{=}\sum _{\mathit{s}}\mathrm{\Delta }{\mathit{T}}_{\mathit{s}}{\mathit{w}}_{\mathit{s}}^{\mathit{T}}{\mathrm{\Omega }}_{\mathit{s}}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}{\mathit{\theta }}_{\mathit{s}}\mathrm{\right)}\mathit{,}$(A.1)where the sum is over the pixels in the map. For all of the methods presented in the paper the weighting per pixel is uniform, except for the apodisation of the different Galactic and point source masks referred to in the main text.

The coupling matrix appearing in Eq. (2) is given by Hivon et al. (2002): $\begin{array}{ccc}{\mathit{M}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{TT}}& \mathrm{=}& \\ & \mathrm{\equiv }& \mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}\begin{array}{c}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{\right)}\end{array}$(A.2)where for the cross spectrum (i,j), $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{W}}_{\mathit{\ell }}\end{array}$ is the power spectrum of the window function $\stackrel{}{{{\mathit{W}}_{\mathit{\ell }}^{\mathit{ij}}}_{\mathrm{˜}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\sum _{\mathit{m}}\stackrel{}{{{\mathit{w}}_{\mathit{\ell m}}^{\mathit{i}}}_{\mathrm{˜}}}\stackrel{}{{{\mathit{w}}_{\mathit{\ell m}}^{\mathit{j}\mathrm{\ast }}}_{\mathrm{˜}}}\mathit{.}$(A.3)

Appendix A.2: Pseudo-C_ℓ covariance matrices

For the case of narrow window functions and uncorrelated pixel noise $\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}$, the covariance matrices can be approximated as $\begin{array}{ccc}\mathrm{⟨}\mathrm{\Delta }\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{T}}_{\mathit{ij}}}}_{\mathrm{˜}}}\mathrm{\Delta }\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{{\mathit{T}}_{\mathit{pq}}}}_{\mathrm{˜}}}\mathrm{⟩}& \mathrm{\approx }& {\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{T}}\mathrm{[}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{\left(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathrm{\left(}\mathit{ip}\mathrm{\right)}\mathrm{\left(}\mathit{jq}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{\right)}\mathrm{+}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathrm{\left(}\mathit{iq}\mathrm{\right)}\mathrm{\left(}\mathit{jp}\mathrm{\right)}}}_{\mathrm{˜}}}{\mathrm{)}}^{\mathrm{]}}\\ & & \mathrm{+}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathit{T}}{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }\mathrm{[}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathrm{2}\mathit{T}\mathrm{\left(}\mathit{ip}\mathrm{\right)}\mathrm{\left(}\mathit{jq}\mathrm{\right)}}}_{\mathrm{˜}}}{\mathrm{)}}^{}\\ & & \\ & & \mathrm{+}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathrm{2}\mathit{T}\mathrm{\left(}\mathit{jp}\mathrm{\right)}\mathrm{\left(}\mathit{iq}\mathrm{\right)}}}_{\mathrm{˜}}}{\mathrm{)}}^{}\\ & & \mathrm{+}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathrm{2}\mathit{T}\mathrm{\left(}\mathit{jq}\mathrm{\right)}\mathrm{\left(}\mathit{ip}\mathrm{\right)}}}_{\mathrm{˜}}}{\mathrm{)}}^{\mathrm{]}}\\ & & \mathrm{+}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathit{TT}\mathrm{\left(}\mathit{ip}\mathrm{\right)}\mathrm{\left(}\mathit{jq}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{+}{\mathrm{\Xi }}_{\mathit{TT}}\mathrm{(}\mathit{\ell ,}{\mathit{\ell }}^{\mathrm{\prime }}\mathit{,}\stackrel{}{{{\mathit{W}}^{\mathit{TT}\mathrm{\left(}\mathit{iq}\mathrm{\right)}\mathrm{\left(}\mathit{jp}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{)}\mathit{,}\end{array}$(A.4)where Ξ is the matrix defined in Eq. (A.2). The window functions are given by:

$\begin{array}{ccc}& & \\ & & \\ & & \end{array}$where $\begin{array}{ccc}& & \stackrel{}{{{\mathit{w}}_{\mathit{\ell m}}^{\mathrm{\left(}\mathit{ij}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{=}\sum _{\mathit{s}}{\mathit{w}}_{\mathit{s}}^{\mathit{i}}{\mathit{w}}_{\mathit{s}}^{\mathit{j}}{\mathrm{\Omega }}_{\mathit{s}}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}{\mathit{\theta }}_{\mathit{s}}\mathrm{\right)}\mathit{,}\\ & & \stackrel{}{{{\mathit{w}}_{\mathit{\ell m}}^{\mathit{T}\mathrm{\left(}\mathit{ij}\mathrm{\right)}}}_{\mathrm{˜}}}\mathrm{=}\sum _{\mathit{s}}\mathrm{\left(}{\mathit{\sigma }}_{\mathit{s}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}{\mathit{w}}_{\mathit{s}}^{\mathit{i}}{\mathit{w}}_{\mathit{s}}^{\mathit{j}}{\mathrm{\Omega }}_{\mathit{s}}^{\mathrm{2}}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\right(}{\mathit{\theta }}_{\mathit{s}}\mathrm{)}\mathit{.}\end{array}$To avoid cumbersome notation, we have omitted indices from the theoretical spectra appearing in Eq. (A.4). In practice, these spectra include unresolved foreground contributions and are smoothed by the appropriate “beam” transfer functions b^ij, which actually describe the overall transfer function from the sky to the maps¹¹. In addition, these covariance matrices are corrected for the pixel window functions p_ℓ (i.e., covariance matrices such as $\mathrm{⟨}\mathrm{\Delta }\stackrel{}{{{\mathit{C}}_{\mathit{\ell }}^{{\mathit{T}}_{\mathit{ij}}}}_{\mathrm{˜}}}\mathrm{\Delta }\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{{\mathit{T}}_{\mathit{pq}}}}_{\mathrm{˜}}}\mathrm{⟩}$ are divided by ${\mathit{p}}_{\mathit{\ell }}^{\mathrm{2}}{\mathit{p}}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathrm{2}}$).

Finally, the Planck maps have correlated pixel noise. We adopt an approximate, heuristic modification of the formulae above to include this effect in the covariance matrices, based on the power spectra of half-ring difference maps. This procedure is detailed in Appendix A.8.

Appendix A.3: Combining intra-frequency cross-spectra

For Planck, the vector containing the power spectra, and its associated covariance matrix, are both large and so require substantial compression to make the computation of a high-ℓ likelihood fast enough for parameter estimation. As described in Appendix A.7, after we correct for the “effective” calibration factors for each individual detector set, the power spectra at each frequency are consistent to extremely high accuracy. Any remaining residuals have a negligible impact on the cosmological analysis. Thus, we combine the cross-spectra from different detectors within a given frequency combination into a single power spectrum. We do not average across frequency combinations since the unresolved foregrounds depend on frequency. Further compression can be accomplished, if desired, only after unresolved foreground parameters have been determined.

We form the linear combination of individual cross-spectra, for each multipole, in the following way: $\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{k}}\\ \mathit{\ell }\end{array}\mathrm{=}\sum _{\mathit{ij}\mathrm{\subset }\mathit{k,}\hspace{0.17em}\mathit{i}\ne \mathit{j}}\frac{{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{T}{\mathit{T}}_{\mathit{ij}}}{\mathit{y}}_{\mathit{i}}{\mathit{y}}_{\mathit{j}}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{ij}}\\ \mathit{\ell }\end{array}}{\mathrm{\left(}{\mathit{b}}_{\mathit{\ell }}^{\mathit{Tij}}{\mathrm{)}}^{\mathrm{2}}\mathrm{\right(}{\mathit{p}}_{\mathit{\ell }}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}}\mathrm{·}$(A.10)Here the index k denotes the particular frequency cross-spectrum combination (e.g., 100 × 100, 143 × 217), the coefficients y_i denote the multiplicative factors for each map, ${\mathit{b}}_{\mathit{\ell }}^{\mathit{Tij}}$ is the (isotropised) beam transfer function for the map combination ij, and p_ℓ is the isotropised pixel window function ¹². The coefficients α_ij are normalized so that $\sum _{\mathit{ij}\mathrm{\subset }\mathit{k,}\hspace{0.17em}\mathit{i}\ne \mathit{j}}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{T}{\mathit{T}}_{\mathit{ij}}}\mathrm{=}\mathrm{1}\mathit{,} {\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{T}{\mathit{T}}_{\mathit{ii}}}\mathrm{=}\mathrm{0.}$(A.11)How can we determine the coefficients α_ij? A near optimal combination, $\mathit{X̂}\begin{array}{c}\mathit{k}\\ \mathit{\ell }\end{array}$, is given by solving $\sum _{\mathit{pq}}\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{pq}\end{array}\mathit{X̂}\begin{array}{c}\mathit{k}\\ \mathit{\ell }\end{array}\mathrm{=}\sum _{\mathit{pq}}\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{pq}\end{array}\mathit{X̂}\begin{array}{c}\mathit{pq}\\ \mathit{\ell }\end{array}\mathit{,}$(A.12)where $\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{pq}\end{array}$ is the block of the inverse covariance matrix appropriate to the spectrum combination k. If the covariance matrix accurately describes the data, the solution of Eq. (A.12) properly accounts for the correlations between the cross-spectra. Solving Eq. (A.12) requires the inversion of a large matrix, so we adopt a simpler solution by weighting each estimate by the diagonal component of the relevant covariance matrix, e.g., ${\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{T}{\mathit{T}}_{\mathit{ij}}}\mathrm{\propto }\mathrm{1}\mathit{/}\mathrm{Cov}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{ij}}\\ \mathit{\ell }\end{array}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{ij}}\\ \mathit{\ell }\end{array}\mathrm{\right)}\mathit{.}$(A.13)This has the effect of assigning each cross-spectrum equal weight in the signal dominated regime and an inverse variance weighting in the noise dominated regime. This is the correct solution in the noise dominated regime. The analysis of intra-frequency residuals presented in Appendix A.7 shows that in the signal dominated regime we see excess variance (with no obvious dependence on the detector/detector set combination) compared to what we expect from instrument noise alone. This excess variance is small compared to the signal and is caused by residual beam errors, consistent with the beam eigenmode amplitudes discussed in Appendix A.7, that are not included in the covariance matrices. This is our justification for assigning roughly equal weight to the spectra in the signal dominated regime.

When we construct a likelihood from the combined estimates we construct the full covariance matrix including cross-correlations between the various spectra. In this matrix, the cross-correlations in the signal dominated regime are dominated by cosmic variance if different masks are used for different frequencies. If identical masks are used for all frequencies, the cross-correlations in the signal-dominated regime are dominated by the cross correlations between the CMB and unresolved foregrounds, which are included in the analytic covariance matrices and act as a regularizing contribution (see Appendix C).

Appendix A.4: Covariance matrix of combined spectra

The estimates of the Planck cross-spectra are linear combinations of the pseudo-C_ℓ estimates, so their covariance matrices are given by e.g., $\mathrm{Cov}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{k}}\\ \mathit{\ell }\end{array}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{{\mathit{k}}^{\mathrm{\prime }}}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{\right)}\mathrm{=}\begin{array}{c}\sum \end{array}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{T}{\mathit{T}}_{\mathit{ij}}}{\mathit{\alpha }}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{T}{\mathit{T}}_{\mathit{pq}}}\mathrm{Cov}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{ij}}\\ \mathit{\ell }\end{array}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{pq}}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{\right)}\mathit{.}$(A.14)Analytic expressions for these covariance matrices have been given in Efstathiou (2004, 2006); Hamimeche & Lewis (2008), and are described in Appendix A.2. The covariance matrices are computed assuming a fixed fiducial theoretical model including an unresolved foreground model for each frequency combination. Typically, the unresolved foregrounds introduce corrections to the covariance matrices of a few percent in the transition region between signal and noise domination. In addition, we compute the fiducial model by applying appropriate beam functions b^T_ij for each detector combination.

As discussed above, the number of coupling matrices required to compute expressions such as Eq. (A.14) scales as ${\mathit{N}}_{\mathrm{map}}^{\mathrm{4}}$ and so becomes prohibitively expensive as the number of cross-spectra becomes large. However, most of these coupling matrices are similar, differing primarily in the amplitude of the noise levels and in minor respects such as a small number of missing pixels. We can therefore adopt the same masks and weightings for groups of cross-spectra and compute coupling matrices only for distinct combinations. This dramatically reduces the computational burden. A similar approach was adopted by Lewis (2008) to analyse the WMAP 5-year temperature maps.

It is also straightforward to calculate covariance matrices for differences between different averages. If we form two spectra averaging over different detector combinations $\begin{array}{ccc}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{a}}\\ \mathit{\ell }\end{array}& \mathrm{=}& \\ \mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{b}}\\ \mathit{\ell }\end{array}& \mathrm{=}& \end{array}$then the covariance matrix of the difference $\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{a}}\\ \mathit{\ell }\end{array}\mathrm{-}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{b}}\\ \mathit{\ell }\end{array}$ is simply $\begin{array}{ccc}\mathrm{Cov}\mathrm{(}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{a}}\\ \mathit{\ell }\end{array}\mathrm{-}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{b}}\\ \mathit{\ell }\end{array}\mathrm{)}& \mathrm{=}& \sum _{\mathit{i}\ne \mathit{j,p}\ne \mathit{q}}\mathrm{(}{\mathit{\alpha }}^{\mathit{ij}}{\mathit{\alpha }}^{\mathit{pq}}\mathrm{+}{\mathit{\beta }}^{\mathit{ij}}{\mathit{\beta }}^{\mathit{pq}}\mathrm{-}{\mathit{\alpha }}^{\mathit{ij}}{\mathit{\beta }}^{\mathit{pq}}\mathrm{-}{\mathit{\beta }}^{\mathit{ij}}{\mathit{\alpha }}^{\mathit{pq}}\mathrm{)}\\ & & \end{array}$(A.17)As a pedagogical illustration of the importance of these correlations, we show in Fig. A.1 two toy power spectra drawn from a Gaussian distribution with mean given by a ΛCDM spectrum and covariance given by $\mathrm{Cov}\mathrm{(}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{a}}\\ \mathit{\ell }\end{array}\mathrm{-}\mathit{Ĉ}\begin{array}{c}{\mathit{T}}_{\mathit{b}}\\ \mathit{\ell }\end{array}\mathrm{)}$. That is, these spectra are not computed from a real sky map, but simply drawn directly from the C_ℓ error model, and therefore exclusively contain correlations modelled by the covariance matrix. The apparent “coherent oscillatory features” seen in the difference spectra (bottom panels) are therefore fully described by the CamSpec covariance matrix, accounting for correlated beam and foreground uncertainties, mask-induced coupling etc. The left panel shows a typical realisation, while the case in the right panel is selected as, visually speaking, one of the most “peculiar” within a relatively small set of simulatoins. When assessing the statistical significance of “unexpected features” in the real CMB spectrum, e.g., similar to those seen in Figs. 15 and 16, it is critical to account for these correlations.

Appendix A.5: The “fiducial Gaussian” approximation

We use a likelihood based on the so-called “fiducial Gaussian” approximation. Here we present our justification for this choice, based on an “expansion in covariance” of the exact likelihood in the exact full-sky, isotropic noise case.

Assuming the CMB, noise, and foregrounds are Gaussian, then the probability distribution for the a_lm coefficients of a collection of maps, given a model, is p = e^{− S}, where $\mathit{S}\mathrm{=}\sum _{\mathit{\ell }}\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{(}\mathrm{tr}\mathrm{[}{{C}\ell }^{-1}{{Ĉ}\mathit{\ell }}^{\mathrm{]}}\mathrm{+}\ln {\left|{{C}}_{\mathit{\ell }}\right|}^{\mathrm{)}}\mathit{,}$(A.18)up to a model-independent normalization. Here $\widehat{{{C}}_{\mathit{\ell }}}$ is the matrix of empirical spectra at a given multipole, and Cℓ are their expectation values for the model in question.

Now, a key point to note is that theoretical power spectra typically differ from each other at each ℓ by less than they differ from the observed $\widehat{{{C}}_{\mathit{\ell }}}$, because of cosmic variance and noise in the latter. So we are justified in expanding Eq. (A.18) about a reasonable fiducial model. Considering a single value of ℓ for simplicity, writing ${C}\mathrm{=}{{C}}_{{f}}\mathrm{+}{\Delta }\mathit{,}$(A.19)we obtain $\begin{array}{ccc}\mathit{S}& \mathrm{=}& {\mathit{S}}_{\mathit{f}}\mathrm{+}\mathrm{(}\mathit{l}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{tr}\mathrm{(}\mathrm{-}{{C}f}^{-1}\mathrm{\Delta }{{{C}}_{{f}}}^{-1}{Ĉ}{l}{{{C}}_{{f}}}^{-1}{\mathrm{\Delta }}^{}\\ & & \mathrm{+}\hspace{0.17em}{{{C}}_{{f}}}^{-1}\mathrm{\Delta }{{{C}}_{{f}}}^{-1}\mathrm{\Delta }{{{C}}_{{f}}}^{-1}{{Ĉ}{l}}^{}\\ & & \mathrm{-}\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{2}}{{{C}}_{{f}}}^{-1}\mathrm{\Delta }{{{C}}_{{f}}}^{-1}\mathrm{\Delta }\mathrm{+}\mathit{...})\end{array}$(A.20)to second order in Δ. We may now complete the square in Δ after extracting a term that is small if the fiducial model is accurate. Up to terms independent of Δ, we have, to second order in Δ: $\begin{array}{ccc}{\mathit{S}}_{\mathrm{2}}& \mathrm{=}& \mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{tr}\mathrm{(}{{C}}_{\mathit{f}}^{-1}\mathrm{\Delta }{{{C}}_{\mathit{f}}}^{-1}\mathrm{\Delta }\mathrm{(}{{C}}_{\mathit{f}}^{-1}{Ĉ}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{)}}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathrm{)}\mathrm{tr}\mathrm{(}\mathrm{(}{C}\mathrm{-}{Ĉ}\mathrm{)}{{C}}_{\mathit{f}}^{-1}\mathrm{(}{C}\mathrm{-}{Ĉ}\mathrm{)}{{{C}}_{\mathit{f}}^{-1}}^{\mathrm{)}}\mathit{.}\end{array}$(A.21)Here we have recombined the perturbation and the fiducial model back together, using Eq. (A.19), in the second term to obtain exactly the “fiducial Gaussian” likelihood. The first term is a correction to the fiducial Gaussian likelihood that is typically small if the fiducial model is accurate.

One can motivate neglecting this term by noticing that in its absence the approximate likelihood is unbiased (as the exact one is). One trades getting second derivatives exactly right in the vicinity of the fiducial model with getting the position, though not the depth, of the minimum right.

Vectorizing the distinct elements of C − Ĉ (following Appendix A of Hamimeche & Lewis 2009), and recognizing the coefficients as the inverse covariance matrix elements of the spectra under the fiducial model, we obtain $\begin{array}{ccc}{\mathit{S}}_{\mathrm{fid}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathit{Ĉ}\mathrm{1}\mathit{T}\mathrm{-}{\mathit{C}}^{\mathrm{1}\mathit{T}}\mathit{,}\mathit{Ĉ}\mathrm{2}\mathit{T}\mathrm{-}{\mathit{C}}^{\mathrm{2}\mathit{T}}\mathit{,}\mathit{...}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}\mathrm{ℳ̂}-1\mathrm{(}\mathit{Ĉ}\mathrm{1}\mathit{T}\mathrm{-}{\mathit{C}}^{\mathrm{1}\mathit{T}}\mathit{,}\mathit{Ĉ}\mathrm{2}\mathit{T}\mathrm{-}{\mathit{C}}^{\mathrm{2}\mathit{T}}\mathit{,}\mathit{...}{\mathrm{)}}^{\mathrm{T}}\mathit{,}\end{array}$(A.22)where is the fiducial covariance matrix of the spectra, and the upper indices run on the different pairs of frequencies.

This suggests an easy generalisation to the coupled cut sky pseudo-spectra, given our calculation of their covariances in Appendix A.2. We now replace the power spectra above with corresponding appropriate averages of detector cross-spectra. With bold face now denoting spectra laid out as vectors, and the grand fiducial covariance, our final action is: $\begin{array}{ccc}{\mathit{S}}_{\mathrm{fid}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{{Ĉ}}^{\mathrm{1}\mathit{T}}\mathrm{-}{{C}}^{\mathrm{1}\mathit{T}}\mathit{,}{{Ĉ}}^{\mathrm{2}\mathit{T}}\mathrm{-}{{C}}^{\mathrm{2}\mathit{T}}\mathit{,}\mathit{...}\mathrm{)}\\ & & \end{array}$(A.23)Another advantage of the “fiducial Gaussian” approximation is that instrumental uncertainties (calibration errors, beam errors, etc.) do not appear in the inverse covariance, but only in the expression of the theoretical spectra C^T in Eq. (A.23) above.

Note that if we fix the foreground model C^Fk for each spectrum k, together with the calibration coefficients and beam parameters, we can minimise the likelihood (Eq. (A.23)) with respect to a “best-fit” primary CMB spectrum. This “best-fit” spectrum is given by the solution of $\sum _{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}{\mathit{\ell }}^{\mathrm{\prime }}}\mathrm{(}\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\end{array}{\mathrm{)}}^{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}}\mathit{Ĉ}\begin{array}{c}\mathrm{CMB}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{=}\sum _{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}{\mathit{\ell }}^{\mathrm{\prime }}}\mathrm{(}\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\end{array}{\mathrm{)}}^{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}}\mathrm{(}{\mathit{c}}^{{\mathit{k}}^{\mathrm{\prime }}}\mathit{Ĉ}\begin{array}{c}{\mathit{k}}^{\mathrm{\prime }}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{-}\mathit{Ĉ}\begin{array}{c}\mathit{F}{\mathit{k}}^{\mathrm{\prime }}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{)}\mathit{,}$(A.24)where the c^k are spectrum effective calibration factors (see Appendix A.7). The covariance matrix of the estimates $\mathit{Ĉ}\begin{array}{c}\mathrm{CMB}\\ \mathit{\ell }\end{array}$ is given by the inverse of the Fisher matrix:$\mathrm{⟨}\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathrm{CMB}\\ \mathit{\ell }\end{array}\hspace{0.17em}\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathrm{CMB}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}\mathrm{=}{\left(\sum _{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}}\mathrm{(}\mathrm{ℳ̂}\begin{array}{c}-1\\ \mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}\end{array}{\mathrm{)}}^{\mathit{k}{\mathit{k}}^{\mathrm{\prime }}}\right)}^{-1}\mathit{.}$(A.25)

Appendix A.6: Uncertainties on individual detector sets beams and calibrations

Let us consider two detectors (or detector sets) X and Y. Neglecting instrumental noise, the cross-spectrum ${\mathit{C}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{obs}}$ is related to the true one, ${\mathit{C}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{sky}}$, through ${\mathit{C}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{obs}}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{sky}}{\mathit{W}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{eff}\mathit{,}\mathrm{true}}\mathit{,}$(A.26)where ${\mathit{W}}_{\mathit{\ell }}^{\mathit{XY,}\mathrm{eff}\mathit{,}\mathrm{true}}$ is the effective beam window function. Note that because of the optical beam non-circularity and the Planck scanning strategy, ${\mathit{W}}_{\mathit{\ell }}^{\mathit{XY}}\ne {{}^{\mathrm{(}}\mathit{W}_{\mathit{\ell }}^{\mathit{XX}}{{\mathit{W}}_{\mathit{\ell }}^{\mathit{YY}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$ when X ≠ Y, while W^XY = W^YX for any X and Y. In the ℓ range of interest, ${\mathit{W}}_{\mathit{\ell }}^{\mathit{XY}}\mathrm{\ge }\mathrm{0}$, so we denote ${\mathit{W}}^{\mathit{XY}}\mathrm{=}{{}^{\mathrm{(}}\mathit{B}^{\mathit{XY}}{}^{\mathrm{)}}}^{\mathrm{2}}$, following the usual prescription for simple (circular) beam models. In what follows, we drop the XY pair superscript except when they are required for clarity.

Our analyses use the best estimated ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{est}}$ of the sky power spectrum, where the measured ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{obs}}$ is corrected by a nominal effective window ${\mathit{W}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{nom}}$: $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{est}}& \mathrm{=}& {\mathit{C}}_{\mathit{\ell }}^{\mathrm{obs}}\mathit{/}{\mathit{W}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathit{,}\\ & \mathrm{=}& {\mathit{C}}_{\mathit{\ell }}^{\mathrm{sky}}{\mathit{W}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{true}}\mathit{/}{\mathit{W}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{nom}}\mathit{,}\\ & \mathrm{=}& {\mathit{C}}_{\mathit{\ell }}^{\mathrm{sky}}{\mathrm{(}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{true}}\mathit{/}{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{nom}}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{.}\end{array}$(A.27)The ratio ${\mathit{B}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{true}}\mathit{/}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{eff}\mathit{,}\mathrm{nom}}$, which determines the uncertainty on the angular power spectrum due to the beam, is estimated using Monte-Carlo simulations of planet transits.

We estimate ${\mathit{B}}_{\mathit{\ell }}^{\mathrm{mean}}$ and ${\mathit{W}}_{\mathit{\ell }}^{\mathrm{mean}}$ from the Monte Carlo simulations as

$\begin{array}{ccc}& & \\ & & \end{array}$and compute the deviations around the mean ${\mathrm{\Delta }}_{\mathit{\ell }}^{\mathit{i}}\mathrm{=}\ln \mathrm{(}{\mathit{B}}_{\mathit{\ell }}^{\mathit{i}}\mathit{/}{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{mean}}}^{\mathrm{)}}\mathit{.}$(A.30)Since the relative dispersion of the simulated ${\mathit{W}}_{\mathit{\ell }}^{\mathit{i}}$ is small (less than 1%), the deviations are well approximated by ${\mathrm{\Delta }}_{\mathit{\ell }}^{\mathit{i}}\mathrm{\simeq }\frac{\mathrm{1}}{\mathrm{2}}\ln \mathrm{(}{\mathit{W}}_{\mathit{\ell }}^{\mathit{i}}\mathit{/}{{\mathit{W}}_{\mathit{\ell }}^{\mathrm{mean}}}^{\mathrm{)}}\mathit{.}$(A.31)The matrix Δ then has n_MC rows and ℓ_max + 1 columns. Its Singular Value Decomposition (SVD) is given by ${\Delta }\mathrm{=}{MDV}T$(A.32)where M is an orthogonal n_MC × n_MC matrix (i.e., M^TM = MM^T = I_nMC), D is a diagonal matrix with n_MC non-negative eigenvalues, and V is a matrix with ℓ_max + 1 rows whose n_MC columns are orthonormal vectors (i.e., V^TV = I_nMC). Here I_nMC is the identity matrix.

The covariance matrix of the beam deviations is defined as $\begin{array}{ccc}{C}& \mathrm{\equiv }& {{\Delta }}^{\mathrm{T}}{\Delta }\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}\mathrm{)}\\ & \mathrm{=}& {V}{{D}}^{\mathrm{2}}{{V}}^{\mathrm{T}}\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,}\end{array}$(A.33)from which we compute the eigenmode matrix ${E}\mathrm{\equiv }{D}{{V}}^{\mathrm{T}}\mathit{/}\mathrm{(}{\mathit{n}}_{\mathrm{MC}}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$(A.34)using the SVD of Δ. Most of the statistical content of Δ or C is limited to the first few modes n_modes with the largest eigenvalues. We therefore keep only the largest n_modes = 5 of the E matrix. The beam uncertainty for a given spectrum is then given by $\begin{array}{ccc}{\mathit{B}}_{\mathit{\ell }}& \mathrm{=}& {\mathit{B}}_{\mathit{\ell }}^{\mathrm{mean}}\exp {\mathrm{(}{{g}}^{\mathrm{T}}{{E}}^{\mathrm{)}}}_{\mathit{\ell }}\\ & \mathrm{=}& {\mathit{B}}_{\mathit{\ell }}^{\mathrm{mean}}\exp \left(\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{modes}}}{\mathit{g}}^{\mathit{k}}{\mathit{E}}_{\mathit{\ell }}^{\mathit{k}}\right)\end{array}$(A.35)where g is a vector of independant Gaussian variates of unit variance with n_modes elements, and E_k(ℓ) is the kth row of E.

This can be generalised to a set of spectra. Taking three pairs of detector sets a = { UV }, b = { XY } and c = { ZT }, one can write $\begin{array}{ccc}\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}\mathrm{=}\mathrm{(}{{M}}^{\mathit{a}}{{M}}^{\mathit{b}}{{M}}^{\mathit{c}}\mathrm{)}\mathit{.}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{.}& & \end{array}$(A.36)and the covariance matrix is given by $\begin{array}{ccc}{{C}}^{\mathit{abc}}& \mathrm{=}& {\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}}^{\mathrm{T}}\mathit{.}\mathrm{(}{{\Delta }}^{\mathit{a}}{{\Delta }}^{\mathit{b}}{{\Delta }}^{\mathit{c}}\mathrm{)}\\ & \mathrm{=}& {\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)}^{\mathrm{T}}\mathit{.}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{.}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{,}\end{array}$(A.37)The beam errors can therefore be correlated (and in fact are strongly so, see Planck Collaboration VII 2014). In the next Appendix, this general covariance matrix is used to derive the beam error eigenmodes of combined spectra for the CamSpec likelihood.

Appendix A.7: Calibration and beam uncertainties for the CamSpec likelihood

Four effective cross-spectra are used in the CamSpec likelihood, each using an individually-prescribed ℓ-range. For each of the effective power spectra, all eligible mask- and beam-deconvolved cross-spectra $\mathit{Ĉ}\begin{array}{c}\mathit{XY}\\ \mathit{l}\end{array}$ are used, weighted according to $\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\mathrm{=}\sum _{\mathit{XY}}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{XY,p}}\widehat{\hspace{0.17em}\mathit{C}}\begin{array}{c}\mathit{XY}\\ \mathit{\ell }\end{array}\mathit{,}$(A.38)with p labelling the effective spectrum. As described in Appendix A.6, uncertainties in the determination of the HFI effective beams are described in terms of beam eigenmodes, E_k, and distributions of the coresponding eigenvalues. To propagate beam errors into the likelihood, we start by using the eigenvalues, along with the ${\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{XY}}$ weights, to construct an appropriate covariance matrix for the effects of beam errors in $\mathrm{⟨}\mathrm{⟨}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}$: $\mathrm{⟨}\mathrm{⟨}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\hspace{0.17em}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}\mathrm{\approx }\mathrm{4}\sum _{\mathit{ij}}\sum _{\mathit{XY}}\sum _{\mathit{ZW}}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{XY,p}}\hspace{0.17em}{\mathit{\alpha }}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{ZW,p}}\hspace{0.17em}{\mathit{E}}_{\mathit{i}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{\mathit{E}}_{\mathit{j}}^{\mathit{XY}}\mathrm{\left(}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{R}}_{\mathit{ij}}^{\mathit{XY,ZW}}\mathit{,}$(A.39)where ${\mathit{R}}_{\mathit{ij}}^{\mathit{XY,ZW}}$ is the correlation between the ith eigenmode of the XY cross-spectrum with the jth eigenmode of the ZW cross-spectrum. The portion of this matrix corresponding to the ℓ-range used in the likelihood is then extracted, and itself singular-value-decomposed. We keep the first n_effmodes (typically five) eigenmodes ${\mathit{E}}_{\mathit{i}}^{\mathit{p}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$, i = 1,...,n_effmodes, orthogonal over the ℓ-range and normalized such that the sum of their outer product directly approximates the covariance: $\mathrm{⟨}\mathrm{⟨}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\hspace{0.17em}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}\mathrm{\approx }\sum _{\mathit{i}}{\mathit{E}}_{\mathit{i}}^{\mathit{p}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{\mathit{E}}_{\mathit{i}}^{\mathit{p}}\mathrm{\left(}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}$(A.40)The eigenmodes are illustrated in Fig. A.2.

Next we calculate a suitable covariance matrix between the eigenmodes. This requires the (non-diagonal in ℓ) inter-effective-spectrum covariance matrices $\mathrm{⟨}\mathrm{⟨}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\hspace{0.17em}\mathit{Ĉ}\begin{array}{c}\mathit{q}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}$, given by: $\mathrm{⟨}\mathrm{⟨}\mathit{Ĉ}\begin{array}{c}\mathit{p}\\ \mathit{\ell }\end{array}\hspace{0.17em}\mathit{Ĉ}\begin{array}{c}\mathit{q}\\ {\mathit{\ell }}^{\mathrm{\prime }}\end{array}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}\mathrm{\approx }\mathrm{4}\sum _{\mathit{ij}}\sum _{\mathit{XY}}\sum _{\mathit{ZW}}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{XY,p}}\hspace{0.17em}{\mathit{\alpha }}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{ZW,q}}\hspace{0.17em}{\mathit{E}}_{\mathit{i}}^{\mathit{XY}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}\hspace{0.17em}{\mathit{E}}_{\mathit{j}}^{\mathit{ZW}}\mathrm{\left(}{\mathit{\ell }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{R}}_{\mathit{ij}}^{\mathit{XY,ZW}}\mathit{.}$(A.41)Given these matrices, we can “stack” the effective spectra to form a data vector X and form a grand beam-covariance matrix ⟨ ⟨ XX^T ⟩ ⟩ _beam. X has length ${\mathit{n}}_{\mathit{X}}\mathrm{=}{\sum }_{\mathit{p}}\mathrm{(}{\mathit{\ell }}_{\max }^{\mathit{p}}\mathrm{-}{\mathit{\ell }}_{\min }^{\mathit{p}}\mathrm{+}\mathrm{1}\mathrm{)}$. Zero-extending each eigenmode at both ends, and arranging these into a matrix, we can form an n_X by n_effmodes·n_eff matrix of extended eigenmodes, E_eff, where n_eff = 4. Now we imagine approximating ⟨ ⟨ XX^T ⟩ ⟩ _beam as a correlated outer product of E_eff, $\mathrm{⟨}\mathrm{⟨}\mathit{X}\hspace{0.17em}{\mathit{X}}^{\mathrm{T}}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}\mathrm{\approx }{\mathit{E}}_{\mathrm{eff}}{\mathit{M}}_{\mathrm{eff}}{\mathit{E}}_{\mathrm{eff}}^{\mathrm{T}}\mathit{.}$(A.42)Requiring that the covariance be chosen to minimise the summed squared-difference between elements on the two sides yields: ${\mathit{M}}_{\mathrm{eff}}\mathrm{=}{\mathit{E}}_{\mathrm{eff}}^{\mathrm{T}}\hspace{0.17em}\mathrm{⟨}\mathrm{⟨}\mathit{X}{\mathit{X}}^{\mathrm{T}}\mathrm{⟩}{\mathrm{⟩}}_{\mathrm{beam}}\hspace{0.17em}{\mathit{E}}_{\mathrm{eff}}\mathit{.}$(A.43)The ${\mathit{E}}_{\mathit{i}}^{\mathit{p}}\mathrm{\left(}\mathit{\ell }\mathrm{\right)}$’s and the associated covariance M_eff are then passed to the likelihood.

Appendix A.8: Noise model of HFI detector sets

For strictly uncorrelated pixel noise $\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}$, and pixel weights w_i, the noise power spectra for the temperature maps are: $\stackrel{}{{{\mathit{N}}^{\mathit{T}}}_{\mathrm{˜}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{\mathit{i}}\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}{\mathit{w}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathit{i}}^{\mathrm{2}}\mathit{,}$(A.44)with contribution to the pseudo-C_ℓ estimates, for uncorrelated noise, of $\stackrel{}{{{\mathit{N}}^{\mathit{T}}}_{\mathrm{˜}}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}\sum _{\mathit{i}}\mathrm{(}{\mathit{\sigma }}_{\mathit{i}}^{\mathit{T}}{\mathrm{)}}^{\mathrm{2}}{\mathit{w}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{\Omega }}_{\mathit{i}}^{\mathrm{2}}\mathit{.}$(A.45)Values for $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{N}}^{\mathit{T}}\end{array}$ are listed in Table A.1. There is a significant dispersion in the noise properties of the two maps at 100 GHz. At 143 and 217 GHz, the PSB maps have significantly lower noise than the SWB maps by a factor of two, as expected.

The noise spectra for the Planck HFI maps are non-white. The following 7-parameter function $\stackrel{}{{{\mathit{N}}_{\mathit{\ell }}^{\mathrm{fit}}}_{\mathrm{˜}}}\mathrm{=}\mathit{A}{\left(\frac{\mathrm{100}}{\mathit{\ell }}\right)}^{\mathit{\alpha }}\mathrm{+}\frac{\mathit{B}\mathrm{(}\mathit{\ell }\mathit{/}\mathrm{1000}{\mathrm{)}}^{\mathit{\beta }}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{(}\mathit{\ell }\mathit{/}{\mathit{\ell }}_{\mathrm{c}}{\mathrm{)}}^{\mathit{\gamma }}{\mathrm{)}}^{\mathit{\delta }}}\mathit{,}$(A.46)is a flexible parameterization that provides accurate fits for all of the HFI channels. The first term on the left models the excess “1/f”-like noise while the second term models the “bell shaped” noise spectrum at high multipoles introduced by time constant deconvolution applied to the time-ordered data, and the low-pass filter designed to remove high-frequency noise due to demodulation. Estimates of the noise spectra can be computed from difference maps constructed from different half-ring surveys¹³. Examples of fits to noise spectra for the 143GHz and 217GHz channels are shown in Fig. A.3. Note that the 100GHz noise spectra are significantly non-white. At 143 and 217GHz, the deviations from white noise are smaller. Since these cross spectra contribute almost all of the weight in the likelihood at high multipoles, the modelling of non-white noise is not a critical factor in forming an accurate likelihood.

We adopt a heuristic approach to fold departures from white noise into the power spectrum covariance estimates. We define a set of noise weight functions, e.g., ${\mathit{\nu }}_{\mathit{\ell }}^{\mathit{m}}\mathrm{=}\frac{\stackrel{}{{{\mathit{N}}_{\mathit{\ell }}^{\mathrm{fit}\mathit{,}\mathrm{m}}}_{\mathrm{˜}}}}{\begin{array}{c}\mathrm{˜}\\ {\mathit{N}}_{\mathit{\ell }}^{\mathit{m}}\end{array}}\mathrm{·}$(A.47)Wherever a σ² term appears in a covariance matrix, we multiply the appropriate coupling matrix by a factor $\mathrm{(}{\mathit{\nu }}_{\mathit{\ell }}^{\mathit{m}}{\mathit{\nu }}_{{\mathit{\ell }}^{\mathrm{\prime }}}^{\mathit{m}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}$(A.48)This heuristic approach can be partially justified by noting that for isotropic Gaussian noise over the full sky, the distribution of Ĉ_ℓ is given by the inverse Wishart distribution: $\mathit{dP}\mathrm{\left(}\mathit{Ĉ}\mathit{\ell }\mathrm{\right|}{C}\ell \mathrm{)}\mathrm{\propto }\mathrm{|}{{W}}_{\mathit{\ell }}{\mathrm{|}}^{\mathrm{-}\left(\frac{\mathrm{2}\mathit{\ell }\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}}{\mathrm{2}}\right)}\exp \left(\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{Tr}{{{W}}_{\mathit{\ell }}}^{-1}{Ĉ}\mathit{\ell }\right)\mathit{,}$(A.49)where ${Ŵ}\ell \mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{\left(}\begin{array}{c}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\mathrm{+}{\mathit{N}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}\mathrm{\right)}\mathit{,}$(A.50)e.g., Percival & Brown (2006). In this special case, our heuristic correction is exact. Further justification of the accuracy of this heuristic approach comes from direct comparisons with numerical simulations incorporating non-white noise (see Sect. 6) and from the accurate agreement of covariance matrices with the ℓ-by-ℓ scatter measured in all of the cross-spectra used to form the likelihood.

Appendix B: Sky masks

Here we show the various masks, which are mentioned in this paper and summarized in Table 2. Figure B.1 shows the set of diffuse Galactic masks, which we obtained by applying a threshold to an ILC-subtracted 353  GHz temperature map. Figure B.2 shows the apodised masks, which we combined to form the set of masks used inthe likelihood analysis.

Appendix C: Chance correlations and inter-frequency consistency tests

Here we explicitly show that, even if the foreground contamination is much smaller than the CMB, chance cross-correlations can produce scatter in the inter-frequency power spectra that is large in the signal dominated regime. To see this, consider the case of two frequencies. Frequency 1 provides a faithful map of the CMB fluctuations. Frequency 2 contains a foreground component F. We therefore write the maps at the two frequencies as: $\begin{array}{ccc}{{X}}_{\mathrm{1}}& \mathrm{=}& \\ {{X}}_{\mathrm{2}}& \mathrm{=}& \end{array}$with spherical transforms $\begin{array}{ccc}{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{1}}& \mathrm{=}& {\mathit{S}}_{\mathit{\ell m}}\mathit{,}\\ {\mathit{a}}_{\mathit{\ell m}}^{\mathrm{2}}& \mathrm{=}& {\mathit{S}}_{\mathit{\ell m}}\mathrm{+}{\mathit{F}}_{\mathit{\ell m}}\mathit{.}\end{array}$The power spectra of the two maps are therefore: $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\sum _{\mathit{m}}{\mathit{S}}_{\mathit{\ell m}}{\mathit{S}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}}\mathit{,}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\sum _{\mathit{m}}\mathrm{(}{\mathit{S}}_{\mathit{\ell m}}\mathrm{+}{\mathit{F}}_{\mathit{\ell m}}\mathrm{)}\mathrm{(}{\mathit{S}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{+}{\mathit{F}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{)}\\ & \mathrm{=}& {\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}}\mathrm{+}\mathrm{2}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}\mathrm{\times }\mathrm{F}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{F}}\mathit{,}\end{array}$and the difference between the power spectra is ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{1}}\mathrm{=}\mathrm{2}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}\mathrm{\times }\mathrm{F}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{F}}\mathit{.}$(C.7)If the CMB is uncorrelated with the foreground, the first term will average to zero over a large number of CMB realisations. But we observe only one realisation of the CMB, and so the cross-term will dominate the inter-frequency residuals even if the foreground contamination is much lower than the CMB (C^F ≪ C^CMB). This is the origin of the excess scatter between the 143 and 217 GHz power spectra at low multipoles shown in Fig. 10.

We construct a specific example of this. The upper map in Fig. C.1 shows an ILC map, estimating the CMB, generated from the 100, 143, 217, and 353 GHz maps. The map in the middle panel shows a “fake” 217 GHz map, i.e. the sum of the ILC map and the 857 GHz map scaled in amplitude to match dust emission at 217 GHz. The real 217 GHz map is shown in the lower panel of Fig. C.1. The “fake” 217 GHz map is evidently quite a good match to the real 217 GHz. By rescaling the 857 GHz map to estimate the dust emission at 143 GHz, we can generate a “fake” 143 GHz map in an analogous way. These “fake” maps each contain two components by construction, and so the inter-frequency residuals from these maps will be dominated by the CMB-foreground cross term in Eq. (C.7).

The 143 and 217 GHz power spectrum difference from these fake maps are compared to the 143 − 217 residuals of the real data in Fig. C.2. The magenta points show the same mask1-mask0 double-difference power spectrum between 217 and 143 GHz as shown in Fig. 3 (with mask0 ≡ G22∪PSA82 and mask1 ≡ G35∪PSA82). The only difference here is that the smoothed dust fit of Eq. (9) has been subtracted from the spectra so that the points scatter around zero. There are advantages to using the double difference because: (a) it is insensitive to calibration differences between frequencies; (b) the contrast between dust emission and other foregrounds (point sources/SZ) is stronger in the area of sky defined by mask1 – mask0 and so the double difference power spectrum should be closer to the results from the fake maps, which use only a dust template ¹⁴. The solid green line shows the double difference power spectrum computed from the fake maps. The amplitude of the scatter from the fake maps and the real data are very similar. In fact, there is almost point-by-point agreement between the results from the real data and the fake maps. This provides compelling evidence that the observed inter-frequency scatter at low multipoles is dominated by the CMB-foreground cross term in Eq. (C.7) rather than some mysterious systematic effect in the data.

The blue points in Fig. C.2 show the difference of the 217 and 143 GHz power spectra for mask1. The scatter at multipoles ≲100 is almost identical to the scatter of the purple points, but increases slightly at higher multipoles. This behaviour is expected and is caused by the additional foreground components (CIB/point sources/SZ), which become comparable in amplitude to Galactic dust at multipoles greater than a few hundred.

We would also expect a strong dependence of the inter-frequency residuals with frequency at low multipoles. Diffuse Galactic emission rises steadily in amplitude from 100 GHz to 217 GHz and hence we would expect the inter-frequency scatter to rise as we go up in frequency. This is what we see in the real data (shown in Fig. C.3). Since diffuse Galactic emission is well approximated by the 857 GHz map at all frequencies, we can predict the scatter seen in this figure by scaling 857 GHz to lower frequencies. For 143 − 217 GHz we observe a scatter of 16(μK)² over the multipole range 50 ≤ ℓ ≤ 500. So, from the 857 GHz scalings to lower frequencies, we predict the scatter given in the Table C.1, which is in excellent agreement with the scatter seen in Fig. C.3.

Appendix D: Validity tests

Appendix D.1: Detailed validity checks

In this appendix we show (Figs. D.1–D.6) the distribution of the cosmological and foreground model parameters for the suite of tests described in Sect. 7. We also show in Figs. D.7 and D.8 the correlation matrix between these parameters and the calibration coefficients of each detector.

Appendix D.2: Cosmological parameters from Planck 70 GHz data

For this analysis we implement the pseudo-C_ℓ method described in Hivon et al. (2002) extended to derive both auto- and cross-power spectra from the 70 GHz maps (see, e.g. Polenta et al. 2005, for a comparison between the two estimators). The noise power spectrum and the covariance matrix are computed using 1000 realistic Planck simulations (FFP6, Planck Collaboration I 2014) of both signal and noise maps. The beam window functions are presented in Planck Collaboration IV (2014), and mode-coupling kernels to correct for incomplete sky coverage are computed from formulae analogous to those in Appendix. A.1.

In Fig. D.9 we show the auto- and cross-power spectra computed from the 70 GHz maps, where cross-spectra are obtained by cross-correlating maps from different pairs of horns (there are three such pairs in total). We use these to construct a likelihood at ℓ > 49 by assuming a Gaussian distribution for the band-powers, and include the covariance matrix estimated from simulations. To estimate cosmological parameters we use this likelihood in combination with the Planck low-ℓ likelihood. We marginalise over a single extragalactic foreground parameter, which is a Poisson term ${\mathit{C}}_{\mathrm{900}}^{\mathit{AS}}$ modelling unresolved residual point sources. Figure D.10 shows the resulting parameters, compared to those from CamSpec. Considering the different ℓ range contributing to the two analysis, the parameter distributions are consistent. In Fig. D.11, we show 70 GHz parameters for three choices of the maximum ℓ considered: ℓ_max = 800, 1000, and 1200. The latter two are consistent. Minor discrepancies are displayed at ℓ_max = 800, which can be explained since there is, in this range, no detection of the point source component. In the same figure, we show results using the Plik likelihood for ℓ_max = 1008, which is consistent with Planck 70 GHz over the same ℓ range.

Appendix D.3: Consistency of the Planck low resolution CMB maps

Here we extend the discussion presented in Sect. 8.3. In Fig. D.12 we show the power spectrum of the residual maps, relative to Commander, for NILC, SEVEM, and SMICA. The maximum discrepancy in the range ℓ ≲ 40 (the multipole where noise begins to become non-negligible) is localised at the quadrupole and is less than 20μK², whereas for the range 3 ≤ ℓ ≲ 40 the differences are of order ≈5μK². Overall, the Planck maps are more in agreement among themselves than with WMAP, except perhaps at the quadrupole. The residual map between WMAP and Planck shows power spectrum residuals from ≈10μK² up to ≈40μK² to ℓ ≲ 40. These figures should be compared to the residual estimated from simulated foreground maps, shown to be ≲10μK² at ℓ ≲ 70 in Planck Collaboration XII (2014).

We complement our results with those obtained from the “FFP6” simulations described in Planck Collaboration XII (2014). They consist of 1000 signal plus noise maps processed through each of the four component separation pipelines.

For each Monte Carlo realisation, we follow the same procedure as in the previous section, i.e. smoothing (FWHM = 329.81′) and binning the maps to N_side = 32. We apply this procedure to both the CMB input maps and the output maps derived by the four component separation algorithms. Again, a Gaussian white noise with a variance of 4μK² is added, and the noise covariance matrix is corrected accordingly. Note that the additional white noise is taken into account not only for numerical regularization (to this extent its amplitude may well be lower), but principally because the output instrumental noise processed through component separation and downgraded to low resolution is far from being white. The additional white noise makes the detailed knowledge of the full noise covariance matrix unimportant.

For each realisation and for each component separation code, we compute the power spectrum of the processed map and of the input signal map, using a common mask. When estimating the power spectrum, the input CMB maps are regularised by adding a negligible Gaussian white noise with 0.1μK² variance. We have checked explicitly that when the white noise is added to the component separated maps, we are able to recover the input power spectrum without bias up to ℓ ~ 60 for all the four component separation methods.

Moreover, we compute the power spectrum of the difference maps (output processed map minus input CMB) for each realisation, in order to evaluate the total amount of residual noise. Note that this is not only given by the added regularization noise, but also from the intrinsic noise, albeit small, that is present in the maps.

Figure D.13 shows the average and the 1σ levels of such noise residuals for each of the component separation solutions. We thus see, that the level of such total noise residuals is well below the difference plotted in Fig. D.12. Therefore, we argue that the existing differences between the codes are due to genuine foreground separation residuals.

Appendix E: Dust cleaning using Planck 353 GHz

The WMAP polarisation products are weighted combinations of Ka, Q, and V band 9 year maps. The WMAP analysis mitigates polarised foreground emission using template fitting for the synchrotron and dust emissions. As a template for synchrotron emission, the WMAP K band channel is used, and for dust a polarisation model is used to create a template map (Page et al. 2007). Here we assess the impact on the WMAP polarisation signal when this dust template is replaced by the Planck 353 GHz map, which provides a more direct tracer of the polarisation. Note that this is only for comparison; in all other analyses, we continue to use the WMAP polarisation products as released by the WMAP team (except using the Planck${\mathit{a}}_{\mathit{\ell m}}^{\mathit{TT}}$ map as discussed in Sect. 8).

Foreground cleaned maps can be written as $\begin{array}{ccc}{{m}}_{\mathrm{clean}}\mathrm{=}{{m}}_{\mathit{i}}\mathrm{-}{\mathit{\alpha }}_{\mathit{i}}{{m}}_{\mathrm{synch}}\mathrm{-}{\mathit{\beta }}_{\mathit{i}}{{m}}_{\mathrm{dust}}& & \end{array}$where m = (Q,U) are linear polarisation Stokes parameter maps, and the index i is for Ka, Q, and V bands. Here m_synch is the WMAP 9-year K band map and for m_dust we use either the WMAP dust template or the Planck 353 GHz maps. For each frequency band, the scaling coefficients α_i and β_i are estimated by minimizing the χ²: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{(}\mathit{\alpha ,}\hspace{0.17em}\mathit{\beta }\mathrm{)}\mathrm{=}{{m}}_{\mathrm{clean}}^{\mathrm{t}}{{C}}^{-1}{{m}}_{\mathrm{clean}}$(E.1)where C is the covariance matrix. Following the WMAP analysis, we do not include the signal contribution in C and we only use the diagonal elements of the noise covariance matrix to estimate χ². Scaling coefficients are computed using two different WMAP masks: the “processing mask” that masks a narrow region in the plane of the Galaxy, and the more conservative “P06” mask used for power spectrum estimation and cosmological analysis.

The template coefficients at each channel are shown in Fig. E.1, estimated using the two different Galactic masks. At Q and V band the estimated coefficients are consistent for the two masks; at Ka band the 353 GHz map gives more stable results than the WMAP dust template. The coefficients using the P06 mask are more uncertain, however, as the residual dust signal outside the mask is low, especially for Ka band. We find that the preferred synchrotron coefficient, α, is slightly lower using the 353 GHz map, and the overall χ², shown in Table E.1, is slightly improved using the Planck dust map.

We now test the effect on cosmological parameters, in particular the optical depth to reionisation, using these two different templates. Using the Planck 353 GHz channel as the dust template, with coefficients estimated using the processing mask, lowers the best fit value of τ by about 1σ (see Fig. E.2). we find τ = 0.075 ± 0.013, compared with τ = 0.089 ± 0.013 using the WMAP dust model. This in turn has the effect of lowering A_s, from 3.088 ± 0.025 to 3.061 ± 0.025, but other ΛCDM parameters are not affected. We note though that using template coefficients estimated outside the P06 Galactic mask, the optical depth using the Planck template is lowered by only 0.5σ compared to the WMAP template, indicating some spatial dependence.We conclude that the impact on cosmological parameters from the choice of dust template is small, but not insignificant.

Still, the low-ℓ polarisation from Planck has not been qualified for cosmological analyses yet, and we only used it to reject the hypothesis of a large effect.

All Tables

All Figures