Planck 2015 results
Free Access
Issue
A&A
Volume 594, October 2016
Planck 2015 results
Article Number A11
Number of page(s) 99
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201526926
Published online 20 September 2016

© ESO, 2016

1. Introduction

This paper presents the angular power spectra of the cosmic microwave background (CMB) and the related likelihood functions, calculated from Planck1 2015 data, which consists of intensity maps from the full mission, along with a subset of the polarization data.

The CMB power spectra contain all of the information available if the CMB is statistically isotropic and distributed as a multivariate Gaussian. For realistic data, these must be augmented with models of instrumental noise, of other instrumental systematic effects, and of contamination from astrophysical foregrounds.

The power spectra are, in turn, uniquely determined by the underlying cosmological model and its parameters. In temperature, the power spectrum has been measured over large fractions of the sky by the Cosmic Background Explorer (COBE; Wright et al. 1996) and the Wilkinson Microwave Anistropy Probe (WMAP; Bennett et al. 2013), and in smaller regions by a host of balloon- and ground-based telescopes (e.g., Netterfield et al. 1997; Hanany et al. 2000; Grainge et al. 2003; Pearson et al. 2003; Tristram et al. 2005b; Jones et al. 2006; Reichardt et al. 2009; Fowler et al. 2010; Das et al. 2011, 2014; Keisler et al. 2011; Story et al. 2013). The Planck 2013 power spectrum and likelihood were discussed in Planck Collaboration XV (2014, hereafter Like13.

The distribution of temperature and polarization on the sky is further affected by gravitational lensing by the inhomogeneous mass distribution along the line of sight between the last scattering surface and the observer. This introduces correlations between large and small scales, which can be estimated by computing the expected contribution of lensing to the 4-point function (i.e., the trispectrum). This can in turn be used to determine the power spectrum of the lensing potential, as is done in Planck Collaboration XV (2016) for this Planck release, and to further constrain the cosmological parameters via a separate likelihood function (Planck Collaboration XIII 2016).

Over the last decade, CMB intensity (temperature) has been augmented by linear polarization data (e.g., Kovac et al. 2002; Kogut et al. 2003; Sievers et al. 2007; Dunkley et al. 2009; Pryke et al. 2009; Araujo et al. 2012; Polarbear Collaboration 2014). Because linear polarization is given by both an amplitude and direction, it can, in turn, be decomposed into two coordinate-independent quantities, each with a different dependence on the cosmology (e.g., Seljak 1997; Kamionkowski et al. 1997; Zaldarriaga & Seljak 1997). One, the so-called E mode, is determined by much the same physics as the intensity, and therefore enables an independent measurement of the background cosmology, as well as a determination of some new parameters (e.g., the reionization optical depth). The other polarization observable, the B mode, is only sourced at early times by gravitational radiation, as produced, for example, during an inflationary epoch. The E and B components are also conventionally taken to be isotropic Gaussian random fields, with only E expected to be correlated with intensity. Thus we expect to be able to measure four independent power spectra, namely the three auto-spectra CTT\hbox{$C_\ell^{TT}$}, CEE\hbox{$C_\ell^{EE}$}, and CBB\hbox{$C_\ell^{BB}$}, along with the cross-spectrum CTE\hbox{$C_\ell^{\textit{TE}}$}.

Estimating these spectra from the likelihood requires cleaned and calibrated maps for all Planck detectors, along with a quantitative description of their noise properties. The required data processing is discussed in Planck Collaboration II (2016), Planck Collaboration III (2016), Planck Collaboration IV (2016), Planck Collaboration V (2016), and Planck Collaboration VIII (2016) for the low-frequency instrument (LFI; 30, 44, and 70 GHz) and Planck Collaboration VII (2016) and Planck Collaboration VIII (2016) for the high-frequency instrument (HFI; 100, 143, 217, 353, 585, and 857 GHz). Although the CMB is brightest over 70–217 GHz, the full range of Planck frequencies is crucial to distinguish between the cosmological component and sources of astrophysical foreground emission, present in even the cleanest regions of sky. We therefore use measurements from those Planck bands dominated by such emission as a template to model the foreground in the bands where the CMB is most significant.

This paper presents the CTT\hbox{$C_\ell^{TT}$}, CEE\hbox{$C_\ell^{EE}$}, and CTE\hbox{$C_\ell^{\textit{TE}}$} spectra, likelihood functions, and basic cosmological parameters from the Planck 2015 release. A complete analysis in the context of an extended ΛCDM cosmology of these and other results from Planck regarding the lensing power spectrum results, as well as constraints from other observations, is given in Planck Collaboration XIII (2016). Wider extensions to the set of models are discussed in other Planck 2015 papers; for example, Planck Collaboration XIV (2016) examines specific models for the dark energy component and extensions to general relativity, and Planck Collaboration XX (2016) discusses inflationary models.

This paper shows that the contribution of high- systematic errors to the polarization spectra are at quite a low level (of the order of a few μK2), therefore enabling an interesting comparison of the polarization-based cosmological results with those derived from CTT\hbox{$C_\ell^{TT}$} alone. We therefore discuss the results for CTE\hbox{$C_\ell^{\textit{TE}}$} and CEE\hbox{$C_\ell^{EE}$} at high multipoles. However, the technical difficulties involved with polarization measurements and subsequent data analysis, along with the inherently lower signal-to-noise ratio (especially for B modes), thus require a careful understanding of the random noise and instrumental and astrophysical systematic effects. For this reason, at large angular scales (i.e., low multipoles ) the baseline results use only a subset of Planck polarization data.

Because of these different sensitivities to systematic errors at different angular scales, as well as the increasingly Gaussian behaviour of the likelihood function at smaller angular scales, we adopt a hybrid approach to the likelihood calculation (Efstathiou 2004, 2006), splitting between a direct calculation of the likelihood on large scales and the use of pseudo-spectral estimates at smaller scales, as we did for the previous release.

The plan of the paper reflects this hybrid approach along with the importance of internal tests and cross-validation. In Sect. 2, we present the low-multipole (< 30) likelihood and its validation. At these large scales, we compute the likelihood function directly in pixel space; the temperature map is obtained by a Gibbs sampling approach in the context of a parameterized foreground model, while the polarized maps are cleaned of foregrounds by a template removal technique.

In Sect. 3, we introduce the high-multipole ( ≥ 30) likelihood and present its main results. At these smaller scales, we employ a pseudo-C approach, beginning with a numerical spherical harmonic transform of the full-sky map, debiased and deconvolved to account for the mask and noise.

Table 1

Likelihood codes and datasets.

Section 4 is devoted to the detailed assessment of this high- likelihood. One technical difference between Like13 and the present work is the move from the CamSpec code to Plik for high- results as well as the released software (Planck Collaboration 2015). The main reason for this change is that the structure of Plik allows more fine-grained tests on the polarization spectra for individual detectors or subsets of detectors. We are able to compare the effect of different cuts on Planck and external data, as well as using methods that take different approaches to estimate the maximum-likelihood spectra from the input maps; these illustrate the small impact of differences in methodology and data preparation, which are difficult to assess otherwise.

We then combine the low- and high- algorithms to form the full Planck likelihood in Sect. 5, assessing there the choice of = 30 for the hybridization scale and establishing the basic cosmological results from Planck 2015 data alone.

Finally, in Sect. 6 we conclude. A series of Appendices discusses sky masks and gives more detail on the individual likelihood codes, both the released version and a series of other codes used to validate the overall methodology.

To help distinguish the many different likelihood codes, which are functions of different parameters and use different input data, Table 1 summarizes the designations used throughout the text.

2. Low-multipole likelihood

At low multipoles, the current Planck release implements a standard joint pixel-based likelihood including both temperature and polarization for multipoles ≤ 29. Throughout this paper, we denote this likelihood “lowTEB”, while “lowP” denotes the polarization part of this likelihood. For temperature, the formalism uses the CMB maps cleaned with Commander (Eriksen et al. 2004, 2008) maps, while for polarization we use the 70 GHz LFI maps and explicitly marginalize over the 30 GHz and 353 GHz maps taken as tracers of synchrotron and dust emission, respectively (see Sect. 2.3), accounting in both cases for the induced noise covariance in the likelihood.

This approach is somewhat different from the Planck 2013 low- likelihood. As described in Like13, this comprised two nearly independent components, covering temperature and polarization information, respectively. The temperature likelihood employed a Blackwell-Rao estimator (Chu et al. 2005) at ≤ 49, averaging over Monte Carlo samples drawn from the exact power spectrum posterior using Commander . For polarization, we had adopted the pixel-based 9-year WMAP polarization likelihood, covering multipoles ≤ 23 (Bennett et al. 2013).

The main advantage of the exact joint approach now employed is mathematical rigour and consistency to higher , while the main disadvantage is a slightly higher computational expense due to the higher pixel resolution required to extend the calculation to = 29 in polarization. However, after implementation of the Sherman-Morrison-Woodbury formula to reduce computational costs (see Appendix B.1), the two approaches perform similarly, both with respect to speed and accuracy, and our choice is primarily a matter of implementational convenience and flexibility, rather than actual results or performance.

2.1. Statistical description and algorithm

We start by reviewing the general CMB likelihood formalism for the analysis of temperature and polarization at low , as described for instance by Tegmark & de Oliveira-Costa (2001), Page et al. (2007), and in Like13. We begin with maps of the three Stokes parameters { T,Q,U } for the observed CMB intensity and linear polarization in some set of HEALPix 2 (Górski et al. 2005) pixels on the sky. In order to use multipoles cut = 29 in the likelihood, we adopt a HEALPix resolution of Nside = 16 which has 3072 pixels (of area 13.6 deg2) per map; this accommodates multipoles up to max = 3Nside−1 = 47, and, considering separate maps of T, Q, and U, corresponds to a maximum of Npix = 3 × 3072 = 9216 pixels in any given calculation, not accounting for any masking.

After component separation, the data vector may be modelled as a sum of cosmological CMB signal and instrumental noise, mX = sX + nX, where s is assumed to be a set of statistically isotropic and Gaussian-distributed random fields on the sky, indexed by pixel or spherical-harmonic indices (ℓm), with X = { T,E,B } selecting the appropriate intensity or polarization component. The signal fields sX have auto- and cross-power spectra CXY\hbox{$C_{\ell}^{XY}$} and a pixel-space covariance matrix S(C)==2maxXYCXYPXY.\begin{equation} \tens{S}(C_{\ell})=\sum_{\ell=2}^{\ell_{\rm max}} \sum_{XY} C_{\ell}^{XY} \tens{P}_{\ell}^{XY}. \end{equation}(1)Here we restrict the spectra to XY = { TT,EE,BB,TE }, with Nside = 16 pixelization, and PXY\hbox{$\tens{P}^{XY}_{\ell}$} is a beam-weighted sum over (associated) Legendre polynomials. For temperature, the explicit expression is (PTT)i,j=2+14πB2P(i·j),\begin{equation} (\tens{P}_{\ell}^{TT})_{i,j} = \frac{2\ell+1}{4\pi} \, B_\ell^2 \, P_\ell (\vec{\hat{n}}_i \cdot \, \vec{\hat{n}}_j), \end{equation}(2)where nˆi\hbox{$\hat{\vec{n}}_i$} is a unit vector pointing towards pixel i, B is the product of the instrumental beam Legendre transform and the HEALPix pixel window, and P is the Legendre polynomial of order ; for corresponding polarization components, see, e.g., Tegmark & de Oliveira-Costa (2001). The instrumental noise is also assumed to be Gaussian distributed, with a covariance matrix N that depends on the Planck detector sensitivity and scanning strategy, and the full data covariance is therefore M = S + N. With these definitions, the full likelihood expression reads (C)=𝒫(m|C)=12π|M|1/2exp(12mTM-1m),\begin{equation} \mathcal{L}(C_{\ell}) = \mathcal{P}(\vec{m}|C_{\ell})=\frac{1}{2\pi|\tens{M}|^{1/2}} \exp\left(-\frac{1}{2}\vec{m}^{\tens{T}}\,\tens{M}^{-1}\vec{m}\right), \label{eq:pbLike} \end{equation}(3)where the conditional probability \hbox{$\mathcal{P}(\vec{m}|C_{\ell})$} defines the likelihood ℒ(C).

The computational cost of this expression is driven by the presence of the matrix inverse and determinant operations, both of which scale computationally as 𝒪(Npix3)\hbox{$\mathcal{O}(N_{\textrm{pix}}^3)$}. For this reason, the direct approach is only computationally feasible at large angular scales, where the number of pixels is low. In practice, we only analyse multipoles below or equal to cut = 29 with this formalism, requiring maps with Nside = 16. Multipoles between cut + 1 and max are fixed to the best-fit ΛCDM spectrum when calculating S. This division between varying and fixed multipoles speeds up the evaluation of Eq. (3) through the Sherman-Morrison-Woodbury formula and the related matrix determinant lemma, as described in Appendix B.1. This results in an order-of-magnitude speed-up compared to the brute-force computation.

2.2. Low- temperature map and mask

Next, we consider the various data inputs that are required to evaluate the likelihood in Eq. (3), and we start our discussion with the temperature component. As in 2013, we employ the Commander algorithm for component separation. This is a Bayesian Monte Carlo method that either samples from or maximizes a global posterior defined by some explicit parametric data model and a set of priors. The data model adopted for the Planck 2015 analysis is described in detail in Planck Collaboration X (2016), and reads sν(θ)=gνi=1NcompFνi(βi,Δν)ai+j=1NtemplateTνjbjν,\begin{equation} \vec{s}_{\nu}(\theta) = g_{\nu} \sum_{i=1}^{N_{\textrm{comp}}} \tens{F}_\nu^i(\beta_i, \Delta_\nu)\,\vec{a}_{i} + \sum_{j=1}^{N_\textrm{template}}{T}^j_{\nu}\vec{b}^\nu_{j}, \end{equation}(4)where θ denotes the full set of unknown parameters determining the signal at frequency ν. The first sum runs over Ncomp independent astrophysical components including the CMB itself; ai is the corresponding amplitude map for each component at some given reference frequency; βi is a general set of spectral parameters for the same component; gν is a multiplicative calibration factor for frequency ν; Δν is a linear correction of the bandpass central frequency; and the function Fνi(βi,Δν)\hbox{$\tens{F}_\nu^i(\beta_i, \Delta_\nu)$} gives the frequency dependence for component i (which can vary pixel-by-pixel and is hence most generally an Npix × Npix matrix). In the second sum, Tνj\hbox{${T}^j_{\nu}$} is one of a set of Ntemplate correction template amplitudes, accounting for known effects such as monopole, dipole, or zodiacal light, with template maps bjν\hbox{$\vec{b}^{\nu}_j$}.

In 2013, only Planck observations between 30 and 353 GHz were employed in the corresponding fit. In the updated analysis, we broaden the frequency range considerably, by including the Planck 545 and 857 GHz channels, the 9-year WMAP observations between 23 and 94 GHz (Bennett et al. 2013), and the Haslam et al. (1982) 408 MHz survey. We can then separate the low-frequency foregrounds into separate synchrotron, free-free, and spinning-dust components, as well as to constrain the thermal dust temperature pixel-by-pixel. In addition, in the updated analysis we employ individual detector and detector-set maps rather than co-added frequency maps, and this gives stronger constraints on both line emission (primarily CO) processes and bandpass measurement uncertainties. For a comprehensive discussion of all these results, we refer the interested reader to Planck Collaboration X (2016).

For the purposes of the present paper, the critical output from this process is the maximum-posterior CMB temperature sky map, shown in the top panel of Fig. 1. This map is natively produced at an angular resolution of 1deg FWHM, determined by the instrumental beams of the WMAP 23 GHz and 408 MHz frequency channels. In addition, the Commander analysis provides a direct goodness-of-fit measure per pixel in the form of the χ2 map shown in Planck Collaboration X (2016, Fig. 22). Thresholding this χ2 map results in a confidence mask that may be used for likelihood analysis, and the corresponding masked region is indicated in the top panel of Fig. 1 by a gray boundary. Both the map and mask are downgraded from their native HEALPix Nside = 256 pixel resolution to Nside = 16 before insertion into the likelihood code, and the map is additionally smoothed to an effective angular resolution of 440′ FWHM.

thumbnail Fig. 1

Top: Commander CMB temperature map derived from the Planck 2015, 9-year WMAP, and 408 MHz Haslam et al. observations, as described in Planck Collaboration X (2016). The gray boundary indicates the 2015 likelihood temperature mask, covering a total of 7% of the sky. The masked area has been filled with a constrained Gaussian realization. Middle: difference between the 2015 and 2013 Commander temperature maps. The masked region indicates the 2013 likelihood mask, removing 13% of the sky. Bottom: comparison of the 2013 (gray) and 2015 (black) temperature likelihood masks.

thumbnail Fig. 2

Top: comparison of the Planck 2013 (blue points) and 2015 (red points) posterior-maximum low- temperature power spectra, as derived with Commander . Error bars indicate asymmetric marginal posterior 68% confidence regions. For reference, we also show the final 9-year WMAP temperature spectrum in light gray points, as presented by Bennett et al. (2013); note that the error bars indicate symmetric Fisher uncertainties in this case. The dashed lines show the best-fit ΛCDM spectra derived from the respective data sets, including high-multipole and polarization information. Middle: difference between the 2015 and 2013 maximum-posterior power spectra (solid black line). The gray shows the same difference after scaling the 2013 spectrum up by 2.4%. Dotted lines indicate the expected ± 1σ confidence region, accounting only for the sky fraction difference. Bottom: reduction in marginal error bars between the 2013 and 2015 temperature spectra; see main text for explicit definition. The dotted line shows the reduction expected from increased sky fraction alone.

The middle panel of Fig. 1 shows the difference between the Planck 2015 and 2013 Commander maximum-posterior maps, where the gray region now corresponds to the 2013 confidence mask. Overall, there are large-scale differences at the 10 μK level at high Galactic latitudes, while at low Galactic latitudes there are a non-negligible number of pixels that saturate the colour scale of ± 25 μK. These differences are well understood. First, the most striking red and blue large-scale features at high latitudes are dominated by destriping errors in our 2013 analysis, due to bandpass mismatch in a few frequency channels effectively behaving as correlated noise during map making. As discussed in section 3 of Planck Collaboration X (2016) and illustrated in Fig. 2 therein, the most significant outliers have been removed from the updated 2015 analysis, and, consequently, the pattern is clearly visible from the difference map in Fig. 1. Second, the differences near the Galactic plane and close to the mask boundary are dominated by negative CO residuals near the Fan region, at Galactic coordinates (l,b) ≈ (110deg,20deg); by negative free-free residuals near the Gum nebula at (l,b) ≈ (260deg,15deg); and by thermal dust residuals along the plane. Such differences are expected because of the wider frequency coverage and improved foreground model in the new fit. In addition, the updated model also includes the thermal Sunyaev-Zeldovich (SZ) effect near the Coma and Virgo clusters in the northern hemisphere, and this may be seen as a roughly circular patch near the Galactic north pole.

Overall, the additional frequency range provided by the WMAP and 408 MHz observations improves the component separation, and combining these data sets makes more sky effectively available for CMB analysis. The bottom panel of Fig. 1 compares the two χ2-based confidence masks. In total, 7% of the sky is removed by the 2015 confidence mask, compared with 13% in the 2013 version.

The top panel in Fig. 2 compares the marginal posterior low- power spectrum, DC( + 1) / (2π), derived from the updated map and mask using the Blackwell-Rao estimator (Chu et al. 2005) with the corresponding 2013 spectrum (Like13). The middle panel shows their difference. The dotted lines indicate the expected variation between the two spectra, σ, accounting only for their different sky fractions3. From this, we can compute χ2==229(D2015D2013σ)2,\begin{equation} \chi^2 = \sum_{\ell=2}^{29} \left(\frac{D_{\ell}^{2015}-D_{\ell}^{2013}}{\sigma_{\ell}}\right)^2, \end{equation}(5)and we find this to be 21.2 for the current data set. With 28 degrees of freedom, and assuming both Gaussianity and statistical independence between multipoles, this corresponds formally to a probability-to-exceed (PTE) of 82%. According to these tests, the observed differences are consistent with random fluctuations due to increased sky fraction alone.

As discussed in Planck Collaboration I (2016), the absolute calibration of the Planck sky maps has been critically reassessed in the new release. The net outcome of this process was an effective recalibration of +1.2% in map domain, or +2.4% in terms of power spectra. The gray line in the middle panel of Fig. 2 shows the same difference as discussed above, but after rescaling the 2013 spectrum up by 2.4%. At the precision offered by these large-scale observations, the difference is small, and either calibration factor is consistent with expectations.

Finally, the bottom panel compares the size of the statistical error bars of the two spectra, in the form of r(σl+σu)|2013(σl+σu)|20151,\begin{equation} r_{\ell} \equiv \frac{\left.\left(\sigma_{\ell}^{\textrm{l}} + \sigma_{\ell}^{\textrm{u}}\right)\right|_{2013}}{\left.\left(\sigma_{\ell}^{\textrm{l}} + \sigma_{\ell}^{\textrm{u}}\right)\right|_{2015}} - 1, \end{equation}(6)where σu\hbox{$\sigma_{\ell}^{\textrm{u}}$} and σl\hbox{$\sigma_{\ell}^{\textrm{l}}$} denote upper and lower asymmetric 68% error bars, respectively. Thus, this quantity measures the decrease in error bars between the 2013 and 2015 spectra, averaged over the upper and lower uncertainties. Averaging over 1000 ideal simulations and multipoles between = 2 and 29, we find that the expected change in the error bar due to sky fraction alone is 7%, in good agreement with the real data. Note that because the net uncertainty of a given multipole is dominated by cosmic variance, its magnitude depends on the actual power spectrum value. Thus, multipoles with a positive power difference between 2015 and 2013 tend to have a smaller uncertainty reduction than points with a negative power difference. Indeed, some multipoles have a negative uncertainty reduction because of this effect.

For detailed discussions and higher-order statistical analyses of the new Commander CMB temperature map, we refer the interested reader to Planck Collaboration X (2016) and Planck Collaboration XVI (2016).

2.3. 70 GHz polarization low-resolution solution

The likelihood in polarization uses only a subset of the full Planck polarization data, chosen to have well-characterized noise properties and negligible contribution from foreground contamination and unaccounted-for systematic errors. Specifically, we use data from the 70 GHz channel of the LFI instrument, for the full mission except for Surveys 2 and 4, which are conservatively removed because they stand as 3σ outliers in survey-based null tests (Planck Collaboration II 2016). While the reason for this behaviour is not completely understood, it is likely related to the fact that these two surveys exhibit the deepest minimum in the dipole modulation amplitude (Planck Collaboration II 2016; Planck Collaboration IV 2016), leading to an increased vulnerability to gain uncertainties and to contamination from diffuse polarized foregrounds.

To account for foreground contamination, the PlanckQ and U 70 GHz maps are cleaned using 30 GHz maps to generate a template for low-frequency foreground contamination, and 353 GHz maps to generate a template for polarized dust emission (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XXX 2016; Planck Collaboration IX 2016). Linear polarization maps are downgraded from high resolution to Nside = 16 employing an inverse-noise-weighted averaging procedure, without applying any smoothing (Planck Collaboration VI 2016).

The final cleaned Q and U maps, shown in Fig. 3, retain a fraction fsky = 0.46 of the sky, masking out the Galactic plane and the “spur regions” to the north and south of the Galactic centre.

thumbnail Fig. 3

Foreground-cleaned, 70 GHz Q (top) and U (bottom) maps used for the low- polarization part of the likelihood. Each of the maps covers 46% of the sky.

At multipoles < 30, we model the likelihood assuming that the maps follow a Gaussian distribution with known covariance, as in Eq. (3). For polarization, however, we use foreground-cleaned maps, explicitly taking into account the induced increase in variance through an effective noise correlation matrix.

To clean the 70 GHz Q and U maps we use a template-fitting procedure. Restricting m to the Q and U maps (i.e., m ≡ [Q,U]) we write m=11αβ(m70αm30βm353),\begin{equation} \vec{m} = \frac{1}{1-\alpha-\beta}\left ( \vec{m}_{70} - \alpha \vec{m}_{30} - \beta \vec{m}_{353} \right)\!, \end{equation}(7)where m70, m30, and m353 are bandpass-corrected versions of the 70, 30, and 353  GHz maps (Planck Collaboration III 2016; Planck Collaboration VII 2016), and α and β are the scaling coefficients for synchrotron and dust emission, respectively. The latter can be estimated by minimizing the quantity χ2=(1αβ)2mTCS+N-1m,\begin{equation} \chi^2 = (1-\alpha-\beta)^2 \vec{m}^\tens{T} \tens{C}_{\tens{S}+\tens{N}}^{-1} \vec{m}, \label{chi2_alphabeta} \end{equation}(8)where CS+N(1αβ)2mmT=(1αβ)2S(C)+N70.\begin{equation} \tens{C}_{\tens{S}+\tens{N}} \equiv (1-\alpha-\beta)^2 \langle \vec{m} \vec{m}^\tens{T} \rangle = (1-\alpha-\beta)^2 \tens{S}(C_\ell)+ \tens{N}_{70} .\label{CSN_alphabeta} \end{equation}(9)Here N70 is the pure polarization part of the 70 GHz noise covariance matrix4 (Planck Collaboration VI 2016), and C is taken as the Planck 2015 fiducial model (Planck Collaboration XIII 2016). We have verified that using the Planck 2013 model has negligible impact on the results described below. Minimization of the quantity in Eq. (8) using the form of the covariance matrix given in Eq. (9) is numerically demanding, since it would require inversion of the covariance matrix at every step of the minimization procedure. However, the signal-to-noise ratio in the 70 GHz maps is relatively low, and we may neglect the dependence on the α and β of the covariance matrix in Eq. (8) using instead: CS+N=S(C)+N70,\begin{equation} \tens{C}_{\tens{S}+\tens{N}} = \tens{S}(C_\ell)+ \tens{N}_{70} ,\label{CSN_alphabeta2} \end{equation}(10)so that the matrix needs to be inverted only once. We have verified for a test case that accounting for the dependence on the scaling parameters in the covariance matrix yields consistent results. We find α = 0.063 and β = 0.0077, with 3σ uncertainties δα ≡ 3σα = 0.025 and δβ ≡ 3σβ = 0.0022. The best-fit values quoted correspond to a polarization mask using 46% of the sky and correspond to spectral indexes (with 2σ errors) nsynch = −3.16 ± 0.40 and ndust = 1.50 ± 0.16, for synchrotron and dust emission respectively (see Planck Collaboration X 2016, for a definition of the foreground spectral indexes). To select the cosmological analysis mask, the following scheme is employed. We scale to 70 GHz both m30 and m353, assuming fiducial spectral indexes nsynch = −3.2 and ndust = 1.6, respectively. In this process, we do not include bandpass correction templates. From either rescaled template we compute the polarized intensities P=Q2+U2\hbox{$P=\sqrt{Q^2+U^2}$} and sum them. We clip the resulting template at equally spaced thresholds to generate a set of 24 masks, with unmasked fractions in the range from 30% to 80% of the sky. Finally, for each mask, we estimate the best-fit scalings and evaluate the probability to exceed, 𝒫(χ2>χ02)\hbox{$\mathcal{P}(\chi^2 > \chi^2_0)$}, where χ02\hbox{$\chi^2_0$} is the value achieved by minimizing Eq. (8). The fsky = 43% processing mask is chosen as the tightest mask (i.e., the one with the greatest fsky) satisfying the requirement \hbox{$\mathcal{P}>5\%$} (see Fig. 4). We use a slightly smaller mask (fsky = 46%) for the cosmological analysis, which is referred to as the R1.50 mask in what follows.

thumbnail Fig. 4

Upper panels: estimated best-fit scaling coefficients for synchrotron (α) and dust (β), for several masks, whose sky fractions are displayed along the bottom horizontal axis (see text). Lower panel: the probability to exceed, 𝒫(χ2>χ02)\hbox{$\mathcal{P}(\chi^2 > \chi^2_0)$}. The red symbols identify the mask from which the final scalings are estimated, but note how the latter are roughly stable over the range of sky fractions. Choosing such a large “processing” mask ensures that the associated errors are conservative.

We define the final polarization noise covariance matrix used in Eq. (3) as N=1(1αβ)2(N70+δα2m30m30T+δβ2m353m353T).\begin{equation} \tens{N} = \frac{1}{(1-\alpha-\beta)^2}\left ( \tens{N}_{70} + \delta^2_\alpha \vec{m}_{30} \vec{m}^\tens{T}_{30} + \delta^2_\beta \vec{m}_{353} \vec{m}^\tens{T}_{353}\right). \label{lowl_covmat} \end{equation}(11)We use 3σ uncertainties, δα and δβ, to define the covariance matrix, conservatively increasing the errors due to foreground estimation. We have verified that the external (column to row) products involving the foreground templates are sub-dominant corrections. We do not include further correction terms arising from the bandpass leakage error budget since they are completely negligible. Intrinsic noise from the templates also proved negligible.

2.4. Low-Planck power spectra and parameters

We use the foreground-cleaned Q and U maps derived in the previous section along with the Commander temperature map to derive angular power spectra. For the polarization part, we use the noise covariance matrix given in Eq. (11), while assuming only 1 μ diagonal regularization noise for temperature. Consistently, a white noise realization of the corresponding variance is added to the Commander map. By adding regularization noise, we ensure that the noise covariance matrix is numerically well conditioned.

For power spectra, we employ the BolPol code (Gruppuso et al. 2009), an implementation of the quadratic maximum likelihood (QML) power spectrum estimator (Tegmark 1997; Tegmark & de Oliveira-Costa 2001). Figure 5 presents all five polarized power spectra. The errors shown in the plot are derived from the Fisher matrix. In the case of EE and TE we plot the Planck 2013 best-fit power spectrum model, which has an optical depth τ = 0.089, as derived from low- WMAP-9 polarization maps, along with the Planck 2015 best model, which has τ = 0.067 as discussed below5. Since the EE power spectral amplitude scales with τ as τ2 (and TE as τ), the 2015 model exhibits a markedly lower reionization bump, which is a better description of Planck data. There is a 2.7σ outlier in the EE spectrum at = 9, not unexpected given the number of low- multipole estimates involved.

thumbnail Fig. 5

Polarized QML spectra from foreground-cleaned maps. Shown are the 2013 Planck best-fit model (τ = 0.089, dot-dashed) and the 2015 model (τ = 0.067, dashed), as well as the 70 GHz noise bias computed from Eq. (11) (blue dotted).

To estimate cosmological parameters, we couple the machinery described in Sect. 2.1 to cosmomc 6 (Lewis & Bridle 2002). We fix all parameters that are not sampled to their Planck 2015 ΛCDM best-fit value (Planck Collaboration XIII 2016) and concentrate on those that have the greatest effect at low : the reionization optical depth τ, the scalar amplitude As, and the tensor-to-scalar ratio r. Results are shown in Table 2 for the combinations (τ,As) and (τ,As,r).

Table 2

Parameters estimated from the low- likelihood.

It is interesting to disentangle the cosmological information provided by low- polarization from that derived from temperature. Low- temperature mainly contains information on the combination Ase− 2τ, at least at multipoles corresponding to angular scales smaller than the scale subtended by the horizon at reionization (which itself depends on τ). The lowest temperature multipoles, however, are directly sensitive to As. On the other hand, large-scale polarization is sensitive to the combination Asτ2. Thus, neither low- temperature nor polarization can separately constrain τ and As. Combining temperature and polarization breaks the degeneracies and puts tighter constraints on these parameters.

In order to disentangle the temperature and polarization contributions to the constraints, we consider four versions of the low-resolution likelihood.

  • 1.

    The standard version described above, which considers the fullset of T, Q, and U maps, along with their covariance matrix, and is sensitive to the TT, TE, EE, and BB spectra.

  • 2.

    A temperature-only version, which considers the temperature map and its regularization noise covariance matrix. It is only sensitive to TT.

  • 3.

    A polarization-only version, considering only the Q and U maps and the QQ, QU, and UU blocks of the covariance matrix. This is sensitive to the EE and BB spectra.

  • 4.

    A mixed temperature-polarization version, which uses the previous polarization-only likelihood but multiplies it by the temperature-only likelihood. This is different from the standard T,Q,U version in that it assumes vanishing temperature-polarization correlations.

The posteriors derived from these four likelihood versions are displayed in Fig. 6. These plots show how temperature and polarization nicely combine to break the degeneracies and provide joint constraints on the two parameters. The degeneracy directions for cases (2) and (3) are as expected from the discussion above; the degeneracy in case (2) flattens for increasing values of τ because for such values the scale corresponding to the horizon at reionization is pulled forward to > 30. By construction, the posterior for case 4 must be equal to the product of the temperature-only (2) and polarization-only (3) posteriors. This is indeed the case at the level of the two-dimensional posterior (see lower right panel of Fig. 6). It is not immediately evident in the one-dimensional distributions because this property does not survive the final marginalization over the non-Gaussian shape of the temperature-only posterior. It is also apparent from Fig. 6 that EE and BB alone do not constrain τ. This is to be expected, and is due to the inverse degeneracy of τ with As, which is almost completely unconstrained without temperature information, and not to the lack of EE signal. By assuming a sharp prior 109Ase− 2τ = 1.88, corresponding to the best estimate obtained when also folding in the high- temperature information (Planck Collaboration XIII 2016), the polarization-only analysis yields τ=0.051-0.020+0.022\hbox{$\tau = 0.051^{+0.022}_{-0.020}$} (red dashed curve in Fig. 6). The latter bound does not differ much from having As constrained by including TT in the analysis, which yields τ=0.054-0.021+0.023\hbox{$\tau = 0.054^{+0.023}_{-0.021}$} (green curves). Finally, the inclusion of non-vanishing temperature-polarization correlations (blue curves) increases the significance of the τ detection at τ = 0.067 ± 0.023. We have also performed a three-parameter fit, considering τ, As, and r for all four likelihood versions described above, finding consistent results.

thumbnail Fig. 6

Likelihoods for parameters from low- data. Panels 1–3: one-dimensional posteriors for log [1010As], τ, and Ase− 2τ for the several sub-blocks of the likelihood, for cases 1 (blue), 2 (black), 3 (red), and 4 (green) – see text for definitions; dashed red is the same as case 3 but imposes a sharp prior 109Ase− 2τ = 1.88. Panel 4: two-dimensional posterior for log [1010As] and τ for the same data combinations; shading indicates the 68% and 95% confidence regions.

2.5. Consistency analysis

Several tests have been carried out to validate the 2015 low- likelihood. Map-based validation and simple spectral tests are discussed extensively in Planck Collaboration IX (2016) for temperature, and in Planck Collaboration II (2016) for Planck 70 GHz polarization. We focus here on tests based on QML and likelihood analyses, respectively employing spectral estimates and cosmological parameters as benchmarks.

We first consider QML spectral estimates C derived using BolPol . To test their consistency, we consider the following quantity: χh2==2max(CCth)M-1(CCth),\begin{equation} \chi_\textrm{h}^2 = \sum_{\ell=2}^{\ell_\textrm{max}} (C_{\ell} - C_{\ell}^\textrm{th}) \, \tens{M}_{\ell \ell^{\prime}}^{-1} \, (C_{\ell} - C_{\ell}^\textrm{th}) , \label{chi2hred} \end{equation}(12)where M=(CCth)(CCth)\hbox{$\tens{M}_{\ell \ell^{\prime}} = \langle (C_{\ell} - C_{\ell}^\textrm{th}) (C_{\ell^{\prime}} - C_{\ell^{\prime}}^\textrm{th}) \rangle $}, Cth\hbox{$C_{\ell}^\textrm{th}$} represents the fiducial Planck 2015 ΛCDM model, and the average is taken over 1000 signal and noise simulations. The latter were generated using the noise covariance matrix given in Eq. (11). We also use the simulations to sample the empirical distribution for χh2\hbox{$\chi_\textrm{h}^2$}, considering both max = 12 (shown in Fig. 7, along with the corresponding values obtained from the data) and max = 30, for each of the six CMB polarized spectra. We report in Table 3 the empirical probability of observing a value of χh2\hbox{$\chi_\textrm{h}^2$} greater than for the data (hereafter, PTE). This test supports the hypothesis that the observed polarized spectra are consistent with Planck’s best-fit cosmological model and the propagated instrumental uncertainties. We verified that the low PTE values obtained for TE are related to the unusually high (but not intrinsically anomalous) estimates 9 ≤ ≤ 11, a range that does not contribute significantly to constraining τ. For spectra involving B, the fiducial model is null, making this, in fact, a null test, probing instrumental characteristics and data processing independent of any cosmological assumptions.

thumbnail Fig. 7

Empirical distribution of χh2\hbox{$\chi_\textrm{h}^2$} derived from 1000 simulations, for the case max = 12 (see text). Vertical bars reindicate the observed values.

Table 3

Empirical probability of observing a value of χh2\hbox{$\chi_h^2$} greater than that calculated from the data.

In order to test the likelihood module, we first perform a 45deg rotation of the reference frame. This leaves the T map unaltered, while sending Q → −U and UQ (and, hence, E → −B and BE). The sub-blocks of the noise covariance matrix are rotated accordingly. We should not be able to detect a τ signal under these circumstances. Results are shown in Fig. 8 for all the full TQU and the TT+EE+BB sub-block likelihoods presented in the previous section. Indeed, rotating polarization reduces only slightly the constraining power in τ for the TT+EE+BB case, suggesting the presence of comparable power in the latter two. On the other hand, τ is not detected at all when rotating the full T,Q,U set, which includes TE and TB. We interpret these results as further evidence that the TE signal is relevant for constraining τ, a result that cannot be reproduced by substituting TB for TE. These findings appear consistent with the visual impression of the low- spectra of Fig. 5. We have also verified that our results stand when r is sampled.

thumbnail Fig. 8

Posterior for τ for both rotated and unrotated likelihoods. The definition and colour convention of the datasets shown are the same as in the previous section (see Fig. 6), while solid and dashed lines distinguish the unrotated and rotated likelihood, respectively.

As a final test of the 2015 Planck low- likelihood, we perform a full end-to-end Monte Carlo validation of its polarization part. For this, we use 1000 signal and noise full focal plane (FFP8) simulated maps (Planck Collaboration XII 2016), whose resolution has been downgraded to Nside = 16 using the same procedure as that applied to the data. We make use of a custom-made simulation set for the Planck 70 GHz channel, which does not include Surveys 2 and 4. For each simulation, we perform the foreground-subtraction procedure described in Sect. 2.3 above, deriving foreground-cleaned maps and covariance matrices, which we use to feed the low- likelihood. As above, we sample only log [1010As] and τ, with all other parameters kept to their Planck best-fit fiducial values. We consider two sets of polarized foreground simulations, with and without the instrumental bandpass mismatch at 30 and 70 GHz. To emphasize the impact of bandpass mismatch, we do not attempt to correct the polarization maps for bandpass leakage. This choice marks a difference from what is done to real data, where the correction is performed (Planck Collaboration II 2016); thus, the simulations that include the bandpass mismatch effect should be considered as a worst-case scenario. This notwithstanding, the impact of bandpass mismatch on estimated parameters is very small, as shown in Fig. 9 and detailed in Table 4. Even without accounting for bandpass mismatch, the bias is at most 1 / 10 of the final 1σ error estimated from real data posteriors. The Monte Carlo analysis also enables us to validate the (Bayesian) confidence intervals estimated by cosmomc on data by comparing their empirical counterparts observed from the simulations. We find excellent agreement (see Table 4).

thumbnail Fig. 9

Empirical distribution of the mean estimated values for log [1010As] (top) and τ (bottom), derived from 1000 FFP8 simulations (see text). For each simulation, we perform a full end-to-end run, including foreground cleaning and parameter estimation. Blue bars refer to simulations that do not include the instrumental bandpass mismatch, while red bars do. The violet bars flag the overlapping area, while the vertical black lines show the input parameters. We note that the (uncorrected) bandpass mismatch effect hardly changes the estimated parameters.

The validation described above only addresses the limited number of instrumental systematic effects that are modelled in the FFP8 simulations, i.e., the bandpass mismatch. Other systematics may in principle affect the measurement of polarization at large angular scales. To address this issue, we have carried out a detailed analysis to quantify the possible impact of LFI-specific instrumental effects in the 70 GHz map (see Planck Collaboration III 2016, for details). Here we just report the main conclusion of that analysis, which estimates the final bias on τ due to all known instrumental systematics to be at most 0.005, i.e., about 0.25σ, well below the final error budget.

Table 4

Statistics for the empirical distribution of estimated cosmological parameters from the FFP8 simulations.

Table 5

Scalings for synchrotron (α) and dust (β) obtained for WMAP, when WMAP K band and Planck 353 GHz data are used as templates.

2.6. Comparison with WMAP-9 polarization cleaned with Planck 353 GHz

In Like13, we attempted to clean the WMAP-9 low resolution maps using a preliminary version of Planck 353 GHz polarization. This resulted in an approximately 1σ shift towards lower values of τ, providing the first evidence based on CMB observations that the WMAP best-fit value for the optical depth may have been biased high. We repeat the analysis here with the 2015 Planck products. We employ the procedure described in Bennett et al. (2013), which is similar to that described above for Planck 2015. However, in contrast to the Planck 70 GHz foreground cleaning, we do not attempt to optimize the foreground mask based on a goodness-of-fit analysis, but stick to the processing and analysis masks made available by the WMAP team. WMAP’s P06 mask is significantly smaller than the 70 GHz mask used in the Planck likelihood, leaving 73.4% of the sky. Specifically we minimize the quadratic form of Eq. (8), separately for the Ka, Q, and V channels from the WMAP-9 release, but using WMAP-9’s own K channel as a synchrotron tracer rather than Planck 30 GHz7. The purpose of the latter choice is to minimize the differences with respect to WMAP’s own analysis. However, unlike the WMAP-9 native likelihood products, which operate at Nside = 8 in polarization, we use Nside = 16 in Q and U, for consistency with the Planck analysis. The scalings we find are consistent with those from WMAP (Bennett et al. 2013) for α in both Ka and Q. However, we find less good agreement for the higher-frequency V channel, where our scaling is roughly 25% lower than that reported in WMAP’s own analysis8. We combine the three cleaned channels in a noise-weighted average to obtain a three-band map and an associated covariance matrix.

We evaluate the consistency of the low-frequency WMAP and Planck 70 GHz low- maps. Restricting the analysis to the intersection of the WMAP P06 and Planck R1.50 masks (fsky = 45.3%), we evaluate half-sum and half-difference Q and U maps. We then compute the quantity χsd2=mTN-1m\hbox{$\chi^2_\mathrm{sd} = \vec{m}^\mathrm{T} \tens{N}^{-1} \vec{m}$} where m is either the half-sum or the half-difference [Q,U] combination and N is the corresponding noise covariance matrix. Assuming that χsd2\hbox{$\chi^2_\mathrm{sd}$} is χ2 distributed with 2786 degrees of freedom we find a PTE(χ2>χsd2)=1.3×10-5\hbox{$\textrm{PTE}(\chi^2>\chi^2_\mathrm{sd}) = 1.3 \times 10^{-5}$} (reduced χ2 = 1.116) for the half-sum, and PTE = 0.84 (reduced χ2 = 0.973) for the half-difference. This strongly suggests that the latter is consistent with the assumed noise, and that the common signal present in the half-sum map is wiped out in the difference.

We also produce noise-weighted sums of the low-frequency WMAP and Planck 70 GHz low-resolution Q and U maps, evaluated in the union of the WMAP P06 and Planck R1.50 masks (fsky = 73.8%). We compute BolPol spectra for the noise-weighted sum and half-difference combinations. These EE, TE, and BB spectra are shown in Fig. 10 and are evaluated in the intersection of the P06 and R1.50 masks. The spectra also support the hypothesis that there is a common signal between the two experiments in the typical multipole range of the reionization bump. In fact, considering multipoles up to max = 12 we find an empirical PTE for the spectra of the half-difference map of 6.8% for EE and 9.5% for TE, derived from the analysis of 10000 simulated noise maps. Under the same hypothesis, but considering the noise-weighted sum, the PTE for EE drops to 0.8%, while that for TE is below the resolution allowed by the simulation set (PTE < 0.1%). The BB spectrum, on the other hand, is compatible with a null signal in both the noise-weighted sum map (PTE = 47.5%) and the half-difference map (PTE = 36.6%).

thumbnail Fig. 10

BolPol spectra for the noise-weighted sum (black) and half-difference (red) WMAP and Planck combinations. The temperature map employed is always the Commander map described in Sect. 2.2 above. The fiducial model shown has τ = 0.065.

We use the Planck and WMAP map combinations to perform parameter estimates from low- data only. We show here results from sampling log [1010As], τ, and the tensor-to-scalar ratio r, with all other parameters kept to the Planck 2015 best fit (the case with r = 0 produces similar results). Figure 11 shows the posterior probability for τ for several Planck and WMAP combinations. They are all consistent, except the Planck and WMAP half-difference case, which yields a null detection for τ – as it should. As above, we always employ the Commander map in temperature. Table 6 gives the mean values for the sampled parameters, and for the derived parameters zre (mean redshift of reionization) and Ase− 2τ. Results from a joint analysis of the WMAP-based low- polarization likelihoods presented here and the Planck high- likelihood are discussed in Sect. 5.7.1.

thumbnail Fig. 11

Posterior probabilities for τ from the WMAP (cleaned with Planck 353 GHz as a dust template) and Planck combinations listed in the legend. Results are presented for the noise-weighted sum both in the union and the intersection of the two analysis masks. The half-difference map is consistent with a null detection, as expected.

Table 6

Selected parameters estimated from the low- likelihood, for Planck, WMAP and their noise-weighted combination.

3. High-multipole likelihood

At high multipoles (> 29), as in Like13, we use a likelihood function based on pseudo-Cs calculated from Planck HFI data, as well as further parameters describing the contribution of foreground astrophysical emission and instrumental effects (e.g., calibration, beams). Aside from the data themselves, the main advances over 2013 include the use of high- polarization information along with more detailed models of foregrounds and instrumental effects.

Section 3.1 introduces the high- statistical description, Sect. 3.2 describes the data we use, Sects. 3.3 and 3.4 describe foreground and instrumental modelling, and Sect. 3.5 describes the covariance matrix between multipoles and spectra. Section 3.6 validates the overall approach on realistic simulations, while Sect. 3.7 addresses the question of the potential impact of low-level instrumental systematics imperfectly corrected by the DPC processing. The reference results generated with the high multipole likelihood are described in Sect. 3.8. A detailed assessment of these results is presented in Sect. 4.

3.1. Statistical description

Assuming a Gaussian distribution for the CMB temperature anisotropies and polarization, all of the statistical information contained in the Planck maps can be compressed into the likelihood of the temperature and polarization auto- and cross-power spectra. In the case of a perfect CMB observation of the full sky (with spatially uniform noise and isotropic beam-smearing), we know the joint distribution of the empirical temperature and polarization power spectra and can build an exact likelihood, which takes the simple form of an inverse Wishart distribution, uncorrelated between multipoles. For a single power spectrum (i.e., ignoring polarization and temperature cross-spectra between detectors) the likelihood for each multipole simplifies to an inverse χ2 distribution with 2 + 1 degrees of freedom. At high enough , the central limit theorem ensures that the shape of the likelihood is very close to that of a Gaussian distributed variable. This remains true for the inverse Wishart generalization to multiple spectra, where, for each , the shape of the joint spectra and cross-spectra likelihood approaches that of a correlated Gaussian (Hamimeche & Lewis 2008; Elsner & Wandelt 2012). In the simple full-sky case, the correlations are easy to compute (Hamimeche & Lewis 2008), and only depend on the theoretical CMB TT, TE, and EE spectra. For small excursions around a fiducial cosmology, as is the case here given the constraining power of the Planck data, one can show that computing the covariance matrix at a fiducial model is sufficient (Hamimeche & Lewis 2008).

The data, however, differ from the idealized case. In particular, foreground astrophysical processes contribute to the temperature and polarization maps. As we see in Sect. 3.3, the main foregrounds in the frequency range we use are emission from dust in our Galaxy, the clustered and Poisson contributions from the cosmic infrared background (CIB), and radio point sources. Depending on the scale and frequency, foreground emission can be a significant contribution to the data, or even exceed the CMB. This is particularly true for dust near the Galactic plane, and for the strongest point sources. We excise the most contaminated regions of the sky (see Sect. 3.2.2). The remaining foreground contamination is taken into account in our model, using the fact that CMB and foregrounds have different emission laws; this enables them to be separated while estimating parameters.

Foregrounds also violate the Gaussian approximation assumed above. The dust distribution, in particular, is clearly non-Gaussian. Following Like13, however, we assume that outside the masked regions we can neglect non-Gaussian features and assume that, as for the CMB, all the relevant statistical information about the foregrounds is encoded in the spatial power spectra. This assumption is verified to be sufficient for our purposes in Sect. 3.6, where we assess the accuracy of the cosmological parameter constraints in realistic Monte Carlo simulations that include data-based (non-Gaussian) foregrounds.

Cutting out the foreground-contaminated regions from our maps biases the empirical power spectrum estimates. We de-bias them using the PolSpice 9 algorithm (Chon et al. 2004) and, following Like13, we take the correlation between multipoles induced by the mask and de-biasing into account when computing our covariance matrix. The masked-sky covariance matrix is computed using the equations in Like13, which are extended to the case of polarization in Appendix C.1.1. Those equations also take into account the inhomogeneous distribution of coloured noise on the sky using a heuristic approach. The approximation of the covariance matrix that can be obtained from those equations is only valid for some specific mask properties, and for high enough multipoles. In particular, as discussed in Appendix C.1.4, correlations induced by point sources cannot be faithfully described in our approximation. Similarly, Monte Carlo simulations have shown that our analytic approximation loses accuracy around = 30. We correct for both of those effects using empirical estimates from Monte Carlo simulations. The computation of the covariance matrix requires knowledge of both the CMB and foreground power spectra, as well as the map characteristics (beams, noise, sky coverage). The CMB and foreground power spectra are obtained iteratively from previous, less accurate versions of the likelihood.

At this stage, we would thus construct our likelihood approximation by compressing all of the individual Planck detector data into mask-corrected (pseudo-) cross-spectra, and build a grand likelihood using these spectra and the corresponding analytical covariance matrix: ln(|C(θ))=12[C(θ)]TC-1[C(θ)]+const.,\begin{equation} -\ln{\cal L}(\vec{\hat C} | \vec{C}(\theta)) = \frac{1}{2} \left[\vec{\hat C} - \vec{C}(\theta)\right]^{\tens{T}} \tens{C}^{-1} \left[\vec{\hat C} - \vec{C}(\theta)\right] + {\rm const.}, \label{eq:basic-likelihood} \end{equation}(13)where Ĉ is the data vector, C(θ) is the model with parameters θ, and C is the covariance matrix. This formalism enables us to separately marginalize over or condition upon different components of the model vector, separately treating cases such as individual frequency-dependent spectra, or temperature and polarization spectra. Obviously, Planck maps at different frequencies have different constraining powers on the underlying CMB, and following Like13 we use this to impose and assess various cuts to keep only the most relevant data.

We therefore consider only the three best CMB Planck channels, i.e., 100 GHz, 143 GHz, and 217 GHz, in the multipole range where they have significant CMB contributions and low enough foreground contamination after masking; we therefore did not directly include the adjacent channels at 70 GHz and 350 GHz in the analysis. In particular, including the 70 GHz data would not bring much at large scales where the results are already cosmic variance limited, and would entail additional complexity in foreground modelling (synchrotron at large scales, additional radio sources excisions at small scales). The cuts in multipole ranges is be described in detail in Sect. 3.2.4. Further, in order to achieve a significant reduction in the covariance matrix size (and computation time), we compress the data vector (and accordingly the covariance matrix), both by co-adding the individual detectors for each frequency and by binning the combined power spectra. We also co-add the two different TE and ET inter-frequency cross-spectra into a single TE spectrum for each pair of frequencies. This compression is lossless in the case without foregrounds. The exact content of the data vector is discussed in Sect. 3.2.

The model vector C(θ) must represent the content of the data vector. It can be written schematically as Cν×νXY|(θ)=MZW,ν×νXY|(θinst)Cν×νZW,sky|(θ)+Nν×νXY|(θinst),\begin{eqnarray} &&\left. C^{XY}_{\nu \times \nu'} \right|_\ell(\theta) = \left. M^{XY}_{ZW,\nu \times \nu'}\right|_\ell (\theta_{\mathrm{inst}})\ \left. C^{ZW,\mathrm{sky}}_{\nu \times \nu'} \right|_\ell(\theta) + \left. N^{XY}_{\nu \times \nu'} \right|_\ell(\theta_{\mathrm{inst}}),\nonumber\\ &&\left. C^{ZW,\mathrm{sky}}_{\nu \times \nu'} \right|_\ell(\theta) = \left. C^{ZW,\mathrm{cmb}} \right|_\ell(\theta) + \left. C^{ZW,\mathrm{fg}}_{\nu \times \nu'} \right|_\ell(\theta), \label{eq:spectrum-model} \end{eqnarray}(14)where Cν×νXY|(θ)\hbox{$\left. C^{XY}_{\nu \times \nu'} \right|_\ell(\theta)$} is the element of the model vector corresponding to the multipole of the XY cross-spectra (X and Y being either T or E) between the pair of frequencies ν and ν. This element of the model originates from the sum of the microwave emission of the sky, i.e., the CMB (CZW,cmb|(θ)) which does not depend of the pair of frequencies (all maps are in units of Kcmb), and foreground (Cν×νZW,fg|(θ)\hbox{$\left. C^{ZW,\mathrm{fg}}_{\nu \times \nu'} \right|_\ell(\theta)$}). Section 3.3 describes the foreground modelling. The mixing matrix MZW,ν×νXY|(θinst)\hbox{$\left. M^{XY}_{ZW,\nu \times \nu'}\right|_\ell (\theta_{\mathrm{inst}})$} accounts for imperfect calibration, imperfect beam correction, and possible leakage between temperature and polarization. It does depend on the pair of frequencies and can depend on the multipole10 when accounting for imperfect beams and leakages. Finally, the noise term Nν×νXY|(θinst)\hbox{$\left. N^{XY}_{\nu \times \nu'} \right|_\ell(\theta_{\mathrm{inst}})$} accounts for the possible correlated noise in the XY cross-spectra for the pair of frequencies ν × ν. Sections 3.2.3 and 3.4 describe our instrument model.

3.2. Data

The data vector Ĉ in the likelihood equation (Eq. (13)) is constructed from concatenated temperature and polarization components, =(TT,EE,TE),\begin{equation} \vec{\hat{C}} = \left( \vec{\hat{C}}^{TT}, \vec{\hat{C}}^{EE}, \vec{\hat{C}}^{\textit{TE}} \right) , \end{equation}(15)which, in turn, comprise the following frequency-averaged spectra: TT=(100×100TT,143×143TT,143×217TT,217×217TT)EE=(100×100EE,100×143EE,100×217EE,143×143EE,143×217EE,217×217EE)TE=(100×100TE,100×143TE,100×217TE,143×143TE,143×217TE,217×217TE).\begin{eqnarray} \vec{\hat{C}}^{TT} \! &=& \left(\vec{\hat{C}}^{TT}_{100 \times 100}, \vec{\hat{C}}^{TT}_{143 \times 143}, \vec{\hat{C}}^{TT}_{143 \times 217}, \vec{\hat{C}}^{TT}_{217 \times 217}\right)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \\ \vec{\hat{C}}^{EE} \! &= &\left(\vec{\hat{C}}^{EE}_{100 \times 100}, \vec{\hat{C}}^{EE}_{100 \times 143}, \vec{\hat{C}}^{EE}_{100 \times 217}, \vec{\hat{C}}^{EE}_{143 \times 143}, \vec{\hat{C}}^{EE}_{143 \times 217}, \vec{\hat{C}}^{EE}_{217 \times 217}\right) \\ \vec{\hat{C}}^{\textit{TE}} \! &=& \left(\vec{\hat{C}}^{\textit{TE}}_{100 \times 100}, \vec{\hat{C}}^{\textit{TE}}_{100 \times 143}, \vec{\hat{C}}^{\textit{TE}}_{100 \times 217}, \vec{\hat{C}}^{\textit{TE}}_{143 \times 143}, \vec{\hat{C}}^{TE }_{143 \times 217}, \vec{\hat{C}}^{\textit{TE}}_{217 \times 217}\right) . \end{eqnarray}The TT data selection is very similar to Like13. We still discard the 100 × 143 and 100 × 217 cross-spectra in their entirety. They contain little extra information about the CMB, as they are strongly correlated with the high S/N maps at 143 and 217 GHz. Including them, in fact, would only give information about the foreground contributions in these cross-spectra, at the expense of a larger covariance matrix with increased condition number. In TE and EE, however, the situation is different since the overall S/N is significantly lower for all spectra, so a foreground model of comparatively low complexity can be used and it is beneficial to retain all the available cross-spectra.

Table 7

Detector sets used to make the maps for this analysis.

We obtain cross power spectra at the frequencies ν × ν using weighted averages of the individual beam-deconvolved, mask-corrected half-mission (HM) map power spectra, XYν×ν|=(i,j)(ν,ν)wi,jXY|×XYi,j|,\begin{equation} \left. \hat{C}_{\nu \times \nu'}^{XY} \right|_{\ell} = \sum_{(i, j) \in (\nu,\nu')} \left. w^{XY}_{i,j}\right|_{\ell} \times \left.\hat{C}^{XY}_{i, j}\right|_{\ell} , \label{eq:hil_cl_freq_avg} \end{equation}(19)where XY ∈ { TT,TE,EE }, and wi,jXY|\hbox{$\left. w^{XY}_{i,j}\right|_{\ell}$} is the multipole-dependent inverse-variance weight for the detector-set map combination (i,j), derived from its covariance matrix (see Sect. 3.5). For XY = TE, we further add the ET power spectra of the same frequency combination to the sum of Eq. (19); i.e., the average includes the correlation of temperature information from detector-set i and polarization information of detector-set j and vice versa.

We construct the Planck high-multipole likelihood solely from the HFI channels at 100, 143, and 217 GHz. These perform best as they have high S/N combined with manageably low foreground contamination. As in Like13, we only employ 70 GHz LFI data for cross-checks (in the high- regime), while the HFI 353 GHz and 545 GHz maps are used to determine the dust model.

3.2.1. Detector combinations

Table 7 summarizes the main characteristics of individual HFI detector sets used in the construction of the likelihood function. As discussed in Sect. 3.1, the likelihood does not use the cross-spectra from individual detector-set maps; instead, we first combine all those contributing at each frequency to form weighted averages. As in 2013, we disregard all auto-power-spectra as the precision required to remove their noise bias is difficult to attain and even small residuals may hamper a robust inference of cosmological parameters (Like13).

In 2015, the additional data available from full-mission observations enables us to construct nearly independent full-sky maps from the first and the second halves of the mission duration. We constructed cross-spectra by cross-correlating the two half-mission maps, ignoring the half-mission auto-spectra at the expense of a very small increase in the uncertainties. This differs from the procedure used in 2013, when we estimated cross-spectra between detectors or detector-sets, and has the advantage of minimizing possible contributions from systematic effects that are correlated in the time domain.

The main motivation for this change from 2013 is that the correlated noise between detectors (at the same or different frequencies) is no longer small enough to be neglected (see Sect. 3.4.4). And while the correction for the “feature” around = 1800, which was (correctly) attributed to residual 4He-JT cooler lines in 2013 (Planck Collaboration VI 2014), has been improved in the 2015 TOI processing pipeline (Planck Collaboration VII 2016), cross-spectra between the two half-mission periods can help to suppress time-dependent systematics, as argued by Spergel et al. (2015). Still, in order to enable further consistency checks, we also build a likelihood based on cross-spectra between full-mission detector-set maps, applying a correction for the effect of correlated noise. The result illustrates that not much sensitivity is lost with half-mission cross-spectra (see the whisker labelled “DS” in Figs. 35, 36, and C.10).

3.2.2. Masks

Temperature and polarization masks are used to discard areas of the sky that are strongly contaminated by foreground emission. The choice of masks is a trade-off between maximizing the sky coverage to minimize sample variance, and the complexity and potentially insufficient accuracy of the foreground model needed in order to deal with regions of stronger foreground emission. The masks combine a Galactic mask, excluding mostly low Galactic-latitude regions, and a point-source mask. We aim to maximize the sky fraction with demonstrably robust results (see Sect. 4.1.2 for such a test).

Temperature masks are obtained by merging the apodized Galactic, CO, and point-source masks described in Appendix A. In polarization, as discussed in Planck Collaboration Int. XXX (2016), even at 100 GHz foregrounds are dominated by the dust emission, so for polarization analysis we employ the same apodized Galactic masks as we use for temperature, because they are also effective in reducing fluctuations in polarized dust emission at the relatively small scales covered by the high- likelihood (contrary to the large Galactic scales), but we do not include a compact-source mask because polarized emission from extragalactic foregrounds is negligible at the frequencies of interest (Naess et al. 2014; Crites et al. 2015).

Table 8

Masks used for the high- analysis.

Table 8 lists the masks used in the likelihood at each frequency channel. We refer throughout to the masks by explicitly indicating the percentage of the sky they retain: T66, T57, T47 for temperature and P70, P50, P41 for polarization. G70, G60, G50, and G41 denote the apodized Galactic masks. As noted above, the apodized P70, P50, and P41 polarization masks are identical to the G70, G50, and G41 Galactic masks.

The Galactic masks are obtained by thresholding the smoothed, CMB-cleaned 353 GHz map at different levels to obtain different sky coverage. All of the Galactic masks are apodized with a 4.̊71 FWHM (σ = 2°) Gaussian window function to localize the mask power in multipole space. In order to adapt to the different relative strengths of signal, noise, and foregrounds, we use different sky coverage for temperature and polarization, ranging in effective sky fraction from 41% to 70% depending on the frequency. The Galactic masks are shown in Fig. 12.

thumbnail Fig. 12

Top: apodized Galactic masks: G41 (blue), G50 (purple), G60 (red), and G70 (orange); these are identical to the polarization masks P41 (used at 217 GHz), P50 (143 GHz), P70 (100 GHz). Bottom: extragalactic-object masks for 217 GHz (purple), 143 GHz (red), and 100 GHz (orange); the CO mask is shown in yellow.

For temperature we use the G70, G60, and G50 Galactic masks at (respectively) 100 GHz, 143 GHz, and 217 GHz. For the first release of Planck cosmological data (Planck Collaboration XI 2016) we made more conservative choices of masks than in this paper (fsky = 49%,31%, and 31% at, respectively, 100, 143, and 217 GHz, to be compared to fsky = 66%,57%, and 47%). Admitting more sky into the analysis requires a thorough assessment of the robustness of the foreground modelling, and in particular of the Galactic dust model (see Sect. 3.3). When retaining more sky close to the Galactic plane at 100 GHz, maps start to show contamination by CO emission that also needs to be masked. This was not the case in the Planck 2013 analysis. We therefore build a CO mask as described in Appendix A. Once we apply this mask, the residual foreground at 100 GHz is consistent with dust and there is no evidence for other anisotropic foreground components, as shown by the double-difference spectra between the 100  GHz band and the 143 GHz band where there is no CO line (Sect. 3.3.1). We also use the CO mask at 217 GHz, although we expect it to have a smaller impact since at this frequency CO emission is fainter and the applied Galactic cut wider. The extragalactic “point” source masks in fact include both point sources and extended objects; they are used only with the temperature maps. Unlike in 2013, we use a different source mask for each frequency, taking into account different source selection and beam sizes (see Appendix A). Both the CO and the extragalactic object masks are apodized with a 30 FWHM Gaussian window function. The different extragalactic masks, as well as the CO mask, are shown in Fig. 12. The resulting mask combinations for temperature are shown in Fig. 13.

thumbnail Fig. 13

Top to bottom: temperature masks for 100 GHz (T66), 143 GHz (T57) and 217 GHz (T47). The colour scheme is the same as in Fig. 12.

3.2.3. Beam and transfer functions

The response to a point source is given by the combination of the optical response of the Planck telescope and feed-horns (the optical beam) with the detector time response and electronic transfer function (whose effects are partially removed during the TOI processing). This response pattern is referred to as the “scanning beam”. It is measured on planet transits (Planck Collaboration VIII 2016). However, the value in any pixel resulting from the map-making operation comes from a sum over many different elements of the timeline, each of which has hit the pixel in a different location and from a different direction. Furthermore, combined maps are weighted sums of individual detectors. All of these result in an “effective beam” window function encoding the multiplicative effect on the angular power spectrum. We note that beam non-circularity and the non-uniform scanning of the sky create differences between auto- and cross-detector beam window functions (Planck Collaboration VII 2014).

In the likelihood analysis, we correct for this by using the effective beam window function corresponding to each specific spectrum; the window functions are calculated with the QuickBeam pipeline, except for one of the alternative analyses (Xfaster ) which relied on the FEBeCoP window functions (see Planck Collaboration VII 2016; Planck Collaboration VII 2014, and references therein for details of these two codes). In Sect. 3.4.3 we discuss the model of their uncertainties.

3.2.4. Multipole range

Following the approach taken in Like13, we use specifically tailored multipole ranges for each frequency-pair spectrum. In general, we exclude multipoles where either the S/N is too low for the data to contribute significant constraints on the CMB, or the level of foreground contamination is so high that the foreground contribution to the power spectra cannot be modelled sufficiently accurately; high foreground contamination would also require us to consider possible non-Gaussian terms in the estimation of the likelihood covariance matrix. We impose the same cuts for the detector-set and half-mission likelihoods for comparison, and we exclude the > 1200 range for the 100 × 100 spectra, where the correlated noise correction is rather uncertain.

thumbnail Fig. 14

Unbinned S/N per frequency for TT (solid blue, for those detector combinations used in the estimate of the TT spectrum), EE (solid red), and TE (solid green). The horizontal orange line corresponds to S/N = 1. The dashed lines indicate the S/N in a cosmic-variance-limited case, obtained by forcing the instrumental noise terms to zero when calculating the power spectrum covariance matrix. The dotted lines indicate the cosmic-variance-limited case computed with the approximate formula of Eq. (20).

thumbnail Fig. 15

Planck power spectra (not yet corrected for foregrounds) and data selection. The coloured tick marks indicate the -range of the cross-spectra included in the Planck likelihood. Although not used in the high- likelihood, the 70 GHz spectra at > 29 illustrate the consistency of the data. The grey line indicates the best-fit Planck 2015 spectrum. The TE and EE plots have a logarithmic horizontal scale for < 30.

Figure 14 shows the unbinned S/N per frequency for TT, EE, and TE, where the signal is given by the frequency-dependent CMB and foreground power spectra, while the noise term contains contributions from cosmic variance and instrumental noise and is given by the diagonal elements of the power-spectrum covariance matrix. The figure also shows the S/N assuming only cosmic variance (CV) in the noise term, obtained either by a full calculation of the covariance matrix with instrumental noise set to zero, or using the approximation σCV{TT,EE}=(2(2+1)fsky)(C{TT,EE})2σCVTE=(2(2+1)fsky)(CTE)2+CTTCEE2·\begin{eqnarray} \label{eq:cv} &&\sigma_{CV}^{\{TT,EE\}} = \sqrt{\left(\frac{2}{(2\ell+1)\fsky}\right)\left(C^{\{TT,EE\}}_\ell\right)^2}\nonumber\\ &&\sigma_{CV}^{\textit{TE}} = \sqrt{\left(\frac{2}{(2\ell+1)\fsky}\right)\frac{\left(C^{\textit{TE}}_\ell\right)^2+C^{TT}_\ell C^{EE}_\ell}{2}} \cdot \end{eqnarray}(20)(see e.g. Percival & Brown 2006).

This figure illustrates that the multipole cuts we apply ensure that the | S/N | ≳ 1. The TT multipole cuts are similar to those adopted in Like13. While otherwise similar to the 2013 likelihood, the revised treatment of dust in the foreground model enables the retention of multipoles < 500 of the 143 × 217 and 217 × 217 GHz TT spectra. As discussed in detail in Sect. 3.3.1, we are now marginalizing over a free amplitude parameter of the dust template, which was held constant for the 2013 release. Furthermore, the greater sky coverage at 100 GHz maximizes its weight at low , so that the best estimate of the CMB signal on large scales is dominated by 100 GHz data. We do not detect noticeable parameter shifts when removing or including multipoles at < 500. See Sect. 4.1 for an in-depth analysis of the impact of different choices of multipole ranges on cosmological parameters.

Table 9

Multipole cuts for the Plik temperature and polarization spectra at high .

thumbnail Fig. 16

The relative weights of each frequency cross-spectrum in the TT (top), TE (middle) and EE (bottom) best-fit solution. Sharp jumps are due to the multipole selection. Weights are normalized to sum to one.

For TE and EE we are more conservative, and cut the low S/N 100 GHz data at small scales (> 1000), and the possibly dust-contaminated 217 GHz at large scales (< 500). Only the 143 × 143TE and EE spectra cover the full multipole range, restricted to < 2000. Retaining more multipoles would require more in-depth modelling of residual systematic effects, which is left to future work. All the cuts are summarized in Table 9 and shown in Fig. 15.

Figure 14 also shows that each of the TT frequency power spectra is cosmic-variance dominated in a wide interval of multipoles. In particular, if we define as cosmic-variance dominated the ranges of multipoles where cosmic variance contributes more than half of the total variance, we find that the 100 × 100 GHz spectrum is cosmic-variance dominated at ≲ 1156, the 143 × 143 GHz at ≲ 1528, the 143 × 217 GHz at ≲ 1607, and the 217 × 217 GHz at ≲ 1566. To determine these ranges, we calculated the ratio of cosmic to total variance, where the cosmic variance is obtained from the diagonal elements of the covariance matrix after setting the instrumental noise to zero. Furthermore, we find that each of the TE frequency power spectra is cosmic-variance limited in some limited ranges of multipoles, below ≲ 150 ( ≲ 50 for the 100 × 100)11, in the range ≈ 250−450 and additionally in the range ≈ 650−700 only for the 100 × 143 GHz and the 143 × 217 GHz power spectra.

Finally, when we co-add the foreground-cleaned frequency spectra to provide the CMB spectra (see Appendix C.4), we find that the CMB TT power spectrum is cosmic-variance dominated at ≲ 1586, while TE is cosmic-variance dominated at ≲ 158 and ≈ 257−464.

Due to the different masks, multipole ranges, noise levels, and to a lesser extent differing foreground contamination, each cross-spectrum ends up contributing differently as a function of scale to the best CMB solution. The determination of the mixing weights is described in Appendix C.4. Figure 16 presents the resulting (relative) weights of each cross-spectra. In temperature, the 100 × 100 spectrum dominates the solution until ≈ 800, when the solution becomes driven by the 143 × 143 up to ≈ 1400. The 143 × 217 and 217 × 217 provide the solution for the higher multipoles. In polarization, the 100 × 143 dominates the solution until ≈ 800 (with an equal contribution from 100 × 100 until ≈ 400 in TE only) while the higher range is dominated by the 143 × 217 contribution. Not surprisingly, the weights of the higher frequencies tend to increase with .

3.2.5. Binning

The 2013 baseline likelihood used unbinned temperature power spectra. For this release, we include polarization, which substantially increases the size of the numerical task. The 2015 likelihood therefore uses binned power spectra by default, downsizing the covariance matrix and speeding up likelihood computations. Indeed, even with the multipole-range cut just described, the unbinned data vector has around 23 000 elements, two thirds of which correspond to TE and EE. For some specific purposes (e.g., searching for oscillatory features in the TT spectrum or testing χ2 statistics) we also produce an unbinned likelihood.

The spectra are binned into bins of width Δ = 5 for 30 ≤ ≤ 99, Δ = 9 for 100 ≤ ≤ 1503, Δ = 17 for 1504 ≤ ≤ 2013, and Δ = 33 for 2014 ≤ ≤ 2508, with a weighting of the C proportional to ( + 1) over the bin widths, Cb==bminbmaxwbC,withwb=(+1)=bminbmax(+1)·\begin{equation} C_{\rm b} = \sum_{\ell=\ell_b^{\mathrm{min}}}^{\ell_b^{\mathrm{max}}} w_b^\ell C_\ell, \ \mathrm{with}\ w_b^\ell = \frac{\ell(\ell+1)}{\sum_{\ell=\ell_b^{\mathrm{min}}}^{\ell_b^{\mathrm{max}}} \ell(\ell+1) } \cdot \end{equation}(21)The bin-widths are odd numbers, since for approximately azimuthal masks we expect a nearly symmetrical correlation function around the central multipole. It is shown explicitly in Sect. 4.1 that the binning does not affect the determination of cosmological parameters in ΛCDM-type models, which have smooth power spectra.

thumbnail Fig. 17

Best foreground model in each of the cross-spectra used for the temperature high- likelihood. The data corrected by the best theoretical CMB C are shown in grey. The bottom panel of each plot shows the residual after foreground correction. The pink line shows the 1σ value from the diagonal of the covariance matrix (32% of the unbinned points are out of this range).

3.3. Foreground modelling

Most of the foreground elements in the model parameter vector are similar to those in Like13. The main differences are in the dust templates, which have changed to accommodate the new masks. The TE and EE foreground model only takes into account the dust contribution and neglects any other Galactic polarized emission, in particular the synchrotron contamination. Nor do we mask out any extragalactic polarized foregrounds, as they have been found to be negligible by ground-based, small-scale experiments (Naess et al. 2014; Crites et al. 2015).

Figure 17 shows the foreground decomposition in temperature for each of the cross-spectra combinations we use in the likelihood. The figure also shows the CMB-corrected data (i.e., data minus the best-fit ΛCDM CMB model) as well as the residuals after foreground correction. In each spectrum, dust dominates the low- modes, while point sources dominate the smallest scales. For 217 × 217 and 143 × 217, the intermediate range has a significant CIB contribution. We note that for 100 × 100, even when including 66% of the sky, the dust contribution is almost negligible and the point-source term is dominant well below = 500. The least foreground-contaminated spectrum is 143 × 143. For comparison, Fig. 18 shows the full model, including the CMB. The foreground contribution is a small fraction of the total power at large scales.

Table 10 summarizes the parameters used for astrophysical foreground modelling and their associated priors.

thumbnail Fig. 18

Best model (CMB and foreground) in each of the cross-spectra used for the temperature high- likelihood. The small light grey points show the unbinned data point, and the dashed grey line show the square root of the noise contribution to the diagonal of the unbinned covariance matrix.

Table 10

Parameters used for astrophysical foregrounds and instrumental modelling.

3.3.1. Galactic dust emission

Galactic dust is the main foreground contribution at large scales and thus deserves close attention. This section describes how we model its power spectra. We express the dust contribution to the power spectrum calculated from map X at frequency ν and map Y at frequency ν as (Cν×νXY,dust)=Aν×νXY,dust×CXY,dust,\begin{equation} \left( C^{XY,\rm dust}_{\nu\times\nu'}\right)_{\ell} = A^{XY,\rm dust}_{\nu\times\nu'}\times C_{\ell}^{XY,\rm dust} , \label{eq:galactic-dust-intensity} \end{equation}(22)where XY is one of TT, EE, or TE, and CXY,dust\hbox{$C_{\ell}^{XY,\rm dust}$} is the template dust power spectrum, with corresponding amplitude Aν×νXY,dust\hbox{$A^{XY,\rm dust}_{\nu\times\nu'}$}. We assume that the dust power spectra have the same spatial dependence across frequencies and masks, so the dependence on sky fraction and frequency is entirely encoded in the amplitude parameter A. We do not try to enforce any a priori scaling with frequency, since using different masks at different frequencies makes determination of this scaling difficult. When both frequency maps ν and ν are used in the likelihood with the same mask, we simply assume that the amplitude parameter can be written as Aν×νXY,dust=aνXY,dust×aνXY,dust.\begin{equation} A^{XY,\rm dust}_{\nu\times\nu'} = a^{XY,\rm dust}_{\nu}\times a^{XY,\rm dust}_{\nu'}. \end{equation}(23)This is clearly not exact when XY = TE and νν. Similarly the multipole-dependent weight used to combine TE and ET for different frequencies breaks the assumption of an invariant dust template. These approximations do not appear to be the limiting factor of the current analysis.

In contrast to the choice we made in 2013, when all Galactic contributions were fixed and a dust template had been explicitly subtracted from the data, we now fit for the amplitude of the dust contribution in each cross-spectrum, in both temperature and polarization. This enables exploration of the possible degeneracy between the dust amplitude and cosmological parameters. A comparison of the two approaches is given in Sect. D.1 and Fig. D.2.

In the following, we describe how we build our template dust power spectrum from high-frequency data and evaluate the amplitude of the dust contamination at each frequency and for each mask.

As we shall see later in Sect. 4.1.2, the cosmological values recovered from TT likelihood explorations do not depend on the dust amplitude priors, as shown by the case “No gal. priors” in Fig. 35 and discussed in Sect. 4.1.2. The polarization case is discussed in Sect. C.3.5. Section 5.3 and Figs. 44 and 45 show the correlation between the dust and the cosmological or other foreground parameters. The dust amplitudes are found to be nearly uncorrelated with the cosmological parameters except for TE. However, the priors do help to break the degeneracies between foreground parameters, which are found to be much more correlated with the dust. In Appendix E we further show that our results are insensitive to broader changes in the dust model.

Galactic TT dust emission.

We use the 545 GHz power spectra as templates for Galactic dust spatial fluctuations. The 353 GHz detectors also have some sensitivity to dust, along with a significant contribution from the CMB, and hence any error in removing the CMB contribution at 353 GHz data translates into biases on our dust template. This is much less of an issue at 545 GHz, to the point where entirely ignoring the CMB contribution does not change our estimate of the template. Furthermore, estimates using 545 GHz maps tend to be more stable over a wider range of multipoles than those obtained from 353 GHz or 857 GHz maps.

We aggressively mask the contribution from point sources in order to minimize their residual, the approximately white spectrum of which is substantially correlated with the value of some cosmological parameters (see the discussion of parameter correlations in Sect. 5.3). The downside of this is that the point-source masks remove some of the brightest Galactic regions that lie in regions not covered by our Galactic masks. This means that we cannot use the well-established power-law modelling advocated in Planck Collaboration XI (2014) and must instead compute an effective dust (residual) template.

All of the masks that we use in this section are combinations of the joint point-source, extended-object, and CO masks used for 100 GHz, 143 GHz, and 217 GHz with Galactic masks of various sizes. In the following discussion we refer only to the Galactic masks, but in all cases the masks contain the other components as well. The half-mission cross-spectra at 545 GHz provide us with a good estimate of the large-scale behaviour of the dust. Small angular scales, however, are sensitive to the CIB, with the intermediate range of scales dominated by the clustered part and the smallest scales by the Poisson distribution of infrared point sources. These last two terms are statistically isotropic, while the dust amplitude depends on the sky fraction. Assuming that the shapes of the dust power spectra outside the masks do not vary substantially as the sky fraction changes, we rely on mask differences to build a CIB-cleaned template of the dust.

Figure 19 shows that this assumption is valid when changing the Galactic mask from G60 to G41. It shows that the 545 GHz cross-half-mission power spectrum can be well represented by the sum of a Galactic template, a CIB contribution, and a point source contribution. The Galactic template is obtained by computing the difference between the spectra obtained in the G60 and the G41 masks. This difference is fit to a simple analytic model CTT,dust(1+hke/t)×(/p)n,\begin{equation} C_{\ell}^{TT,\rm dust} \propto (1 + h \, \ell^k \, {\rm e}^{-\ell/t}) \times (\ell/\ell_p)^n, \label{eq:dust:TTtemplate} \end{equation}(24)with h = 2.3 × 10-11, k = 5.05, t = 56, n = −2.63, and fixing p = 200. The model behaves like a Cℓ,dustTT-2.63\hbox{$C_{\ell,\rm dust}^{\TT} \propto \ell^{-2.63}$} power law at small scales, and has a bump around = 200. The CIB model we use is described in Sect. 3.3.2.

We can compare this template model with the dust content in each of the power spectra we use for the likelihood. Of course those power spectra are strongly dominated by the CMB, so, to reveal the dust content, one has to rely on the same trick that was used for 545 GHz. This however is not enough, since the CMB cosmic variance itself is significant compared to the dust contamination. We can build an estimate of the CMB cosmic variance by assuming that at 100 GHz the dust contamination is small enough that a mask difference gives us a good variance estimate.

thumbnail Fig. 19

Dust model at 545  GHz. The dust template is based on the G60–G41 mask difference of the 545 GHz half-mission cross-spectrum (blue line and circles, rescaled to the dust level in mask G60). Coloured diamonds display the difference between this model (rescaled in each case) and the cross half-mission spectra in the G41, G50, and G60 masks. The residuals are all in good agreement (less so at low , because of sample variance) and are well described by the CIB+point source prediction (orange line). Individual CIB and point sources contributions are shown as dashed and dotted orange lines. The red line is the sum of the dust model, CIB, and point sources for the G60 mask, and is in excellent agreement with the 545 GHz cross half-mission spectrum in G60 (red squares). In all cases, the spectra were computed by using different Galactic masks supplemented by the single combination of the 100 GHz, 143 GHz, and 217 GHz point sources, extended objects and CO masks.

thumbnail Fig. 20

Dust model versus data. In blue, the power spectrum of the double mask difference between 217 GHz and 100 GHz half-mission cross-spectra in masks G60 and G41 (complemented by the joint masks for CO, extended objects, and point sources). In orange, the equivalent spectrum for 143 and 100 GHz. The mask difference enables us to remove the contribution from all the isotropic components (CMB, CIB, and point sources) in the mean. But simple mask differences are still affected by the difference of the CMB in the two masks due to cosmic variance. Removing the 100 GHz mask difference, which is dominated by the CMB, reduces the scatter significantly. The error bars are computed as the scatter in bins of size Δ = 50. The dust model (green) based on the 545 GHz data has been rescaled to the expected dust contamination in the 217 GHz mask difference using values from Table 11. The 143 GHz double mask difference is also rescaled to the level of the 217 GHz difference; i.e., it is multiplied by approximately 14. Different multipole bins are used for the 217 GHz and 143 GHz data to improve readability.

Table 11

Contamination level in each frequency, D = 200.

Figure 20 shows the mask difference (corrected for cosmic variance) between G60 and G41 for the 217 GHz and 143 GHz half-mission cross-spectra, as well as the dust model from Eq. (24). The dust model has been rescaled to the expected mask difference dust residual for the 217 GHz. The 143 GHz mask-difference has also been rescaled in a similar way. The ratio between the two is about 14. Rescaling factors are obtained from Table 11. Error bars are estimated based on the scatter in each bin. The agreement with the model is very good at 217 GHz, but less good at 143 GHz where the greater scatter is probably dominated at large scales by the chance correlation between CMB and dust (which, as we see in Eq. (25), varies as the square root of the dust contribution to the spectra), and at small scale by noise. We also tested these double differences for other masks, namely G50G41 and G60G50, and verified that the results are similar (i.e., general agreement although with substantial scatter).

Finally, we can estimate the level of the dust contamination in each of our frequency maps used for CMB analysis by computing their cross-spectra with the 545 GHz half-mission maps. Assuming that all our maps mν have in common only the CMB and a variable amount of dust, and assuming that m545 = mcmb + a545mdust, the cross-spectra between each of our CMB frequencies maps and the 545 GHz map is (C545×νTT)=CTT,cmb+a545TT,dustaνTT,dustCTT,dust\begin{eqnarray} \left ( C^{TT}_{545\times \nu}\right)_{\ell} &=& C^{TT,\rm cmb}_{\ell} + a^{TT,\rm dust}_{545}a^{TT,\rm dust}_{\nu}\, C^{TT,\rm dust}_{\ell}\nonumber\\ &&\quad + (a^{TT,\rm dust}_{545}+a^{TT,\rm dust}_{\nu})\, C^{\rm chance}_{\ell},\label{eq:hil:chancedust} \end{eqnarray}(25)where Cchance\hbox{$C^{\rm chance}_{\ell}$} is the chance correlation between the CMB and dust distribution (which would vanish on average over many sky realizations). By using the 100 GHz spectrum as our CMB estimate and assuming that the chance correlation is small enough, one can measure the amount of dust in each frequency map by fitting the rescaling factor between the (CMB cleaned) 545 GHz spectrum and the cross frequency spectra. This approach is limited by the presence of CIB which has a slightly different emission law than the dust. We thus limit our fits to the multipoles < 1000 where the CIB is small compared to the dust and we ignore the emission-law differences.

Table 11 reports the results of those fits at each frequency, for each Galactic mask. The error range quoted corresponds to the error of the fits, taking into account the variations when changing the multipole range of the fit from 30 ≤ ≤ 1000 to 30 ≤ ≤ 500. The values reported correspond to the sum of the CIB and the dust contamination at = 200. The last column gives the estimate of the CIB contamination at the same multipole from the joint cosmology and foreground fit. From this table, the ratio of the dust contamination at map level between the 217 GHz and 100 GHz is around 7, while the ratio between the 217 GHz and 143 GHz is close to 3.7.

We derive our priors on the foreground amplitudes from this table, combining the 545 GHz fit with the estimated residual CIB contamination, to obtain the following values: (7 ± 2) μ for the 100 × 100 spectrum (G70); (9 ± 2) μ for 143 × 143 (G60); and (80 ± 20) μ for 217 × 217 (G50). Finally the 143 × 217 value is obtained by computing the geometrical average between the two auto spectra under the worst mask (G60), yielding (21 ± 8.5) μ.

Galactic TE and EE dust emission.

We evaluate the dust contribution in the TE and EE power spectra using the same method as for the temperature. However, instead of the 545 GHz data we use the maps at 353 GHz, our highest frequency with polarization information. At sufficently high sky fractions, the 353 GHz TE and EE power spectra are dominated by dust. As estimated in Planck Collaboration Int. XXX (2016), there is no other significant contribution from the Galaxy, even at 100 GHz. Following Planck Collaboration Int. XXX (2016), and since we do not mask any “point-source-like” region of strong emission, we can use a power-law model as a template for the polarized Galactic dust contribution. Enforcing a single power law for TE and EE and our different masks, we obtain an index of n = −2.4. We use the same cross-spectra-based method to estimate the dust contamination. The dust contribution being smaller in polarization, removing the CMB from the 353 × 353 and the 353 × ν (with ν being one of 100, 143 or 217) is particularly important. Our two best CMB estimates in EE and TE being 100 and the 143 GHz, we checked that using any of 100 × 100, 143 × 143, or 100 × 143 does not change the estimates significantly. Table 12 gives the resulting values. As for the TT case, the cross-frequency, cross-masks estimates are obtained by computing the geometric average of the auto-frequency contaminations under the smallest mask.

Table 12

TE and EE dust contamination levels, D = 500.

3.3.2. Extragalactic foregrounds

The extragalactic foreground model is similar to that of 2013 and in the following we describe the differences. Since we are neglecting any possible contribution in polarization from extragalactic foregrounds, we omit the TT index in the following descriptions of the foreground models. The amplitudes are expressed as \hbox{${\cal D}_\ell$} at = 3000 so that, for any component, the template, C3000FG\hbox{$C^{\rm FG}_{3000}$}, satisfies C3000FG𝒜3000=1\hbox{$C^{\rm FG}_{3000}\, {\cal A}_{3000}=1$} with \hbox{${\cal A}_\ell = \ell(\ell+1)/(2\pi)$}.

The cosmic infrared background.

The CIB model has a number of differences from that used in Like13. First of all, it is now entirely parameterized by a single amplitude 𝒟217CIB\hbox{${\cal D}_{217}^{\rm CIB}$} and a template CCIB\hbox{$C_\ell^{\rm CIB}$}: (Cν×νCIB)=aνCIBaνCIBCCIB×𝒟217CIB,\begin{equation} \left( C^{\rm CIB}_{\nu\times\nu'}\right)_{\ell} = a^{\rm CIB}_{\nu} a^{\rm CIB}_{\nu'} C_{\ell}^{\rm CIB} \times {\cal D}_{217}^{\rm CIB}, \label{eq:CIB} \end{equation}(26)where the spectral coefficients aνCIB\hbox{$a^{\rm CIB}_{\nu}$} represent the CIB emission law normalized at ν = 217  GHz.

In 2013, the template was an effective power-law model with a variable index with expected value n = −1.37 (when including the “highL” data from ACT and SPT). We did not assume any emission law and fitted the 143 GHz and 217 GHz amplitude, along with their correlation coefficient. The Planck Collaboration has studied the CIB in detail in Planck Collaboration XXX (2014) and now proposes a one-plus-two-halo model, which provides an accurate description of the Planck and IRAS CIB spectra from 3000 GHz down to 217 GHz. We extrapolate this model here, assuming it remains appropriate in describing the 143 GHz and 100 GHz data. The CIB emission law and template are computed following Planck Collaboration XXX (2014). The template power spectrum provided by this work has a very small frequency dependence that we ignore.

At small scales, > 2500, the slope of the template is similar to the power law used in Like13. At larger scales, however, the slope is much shallower. This is in line with the variation we observed in 2013 on the power-law index of our simple CIB model when changing the maximum multipole. The current template is shown as the green line in the TT foreground component plots in Fig. 17.

In 2013, the correlation between the 143 GHz and 217 GHz CIB spectra was fitted, favouring a high correlation, greater than 90% (when including the “highL” data). The present model yields a fully correlated CIB between 143 GHz and 217 GHz.

We now include the the CIB contribution at 100 GHz, which was ignored in 2013. Another difference with the 2013 model is that the parameter controlling the amplitude at 217 GHz now directly gives the amplitude in the actual 217 GHz Planck band at = 3000, i.e., it includes the colour correction. The ratio between the two is 1.33. The 2013 amplitude of the CIB contribution at = 3000 (including the highL data) was 66 ± 6.7 μK2, while our best estimate for the present analysis is 63.9 ± 6.6 μK2 (PlanckTT + lowP).

Point sources.

At the likelihood level, we cannot differentiate between the radio- and IR-point sources. We thus describe their combined contribution by their total emissivity per frequency pair, (Cν×νPS)=𝒟ν×νPS/𝒜3000,\begin{equation} \left( C^{\rm PS}_{\nu\times\nu'}\right)_{\ell} = {\cal D}_{\nu\times\nu'}^{\rm PS}/{\cal A}_{3000}, \label{eq:point-sources} \end{equation}(27)where \hbox{${\cal D}_{\nu\times\nu'}$} is the amplitude of the point-source contribution in \hbox{$\mathcal{D}_\ell$} at = 3000. Contrary to 2013, we do not use a correlation parameter to represent the 143 × 217 point-source contribution; instead we use a free amplitude parameter. This has the disadvantage of not preventing a possible unphysical solution. However, it simplifies the parameter optimization, and it is easier to understand in terms of contamination amplitude.

Kinetic SZ (kSZ).

We use the same model as in 2013. The kSZ emission is parameterized with a single amplitude and a fixed template from Trac et al. (2011), (Cν×νkSZ)=CkSZ×𝒟kSZ,\begin{equation} \left( C^{\rm kSZ}_{\nu\times\nu'}\right)_{\ell} = C_{\ell}^{\rm kSZ} \times {\cal D}^\mathrm{kSZ} , \label{eq:kSZ} \end{equation}(28)where \hbox{${\cal D}^\mathrm{kSZ}$} is the kSZ contribution at = 3000.

Thermal SZ (tSZ).

Here again, we use the same model as in 2013. The tSZ emission is also parameterized by a single amplitude and a fixed template using the ϵ = 0.5 model from Efstathiou & Migliaccio (2012), (Cν×νtSZ)=aνtSZaνtSZCtSZ×𝒟143tSZ,\begin{equation} \left( C^{\rm tSZ}_{\nu\times\nu'}\right)_{\ell} = a^{\rm tSZ}_{\nu} a^{\rm tSZ}_{\nu'} C_{\ell}^{\rm tSZ} \times {\cal D}_{143}^\mathrm{tSZ} , \label{eq:tSZ} \end{equation}(29)where aνtSZ\hbox{$ a^{\rm tSZ}_{\nu} $} is the thermal Sunyaev-Zeldovich spectrum, normalized to ν0 = 143  GHz and corrected for the Planck bandpass colour corrections. Ignoring the bandpass correction, we recall that the tSZ spectrum is given by aνtSZ=f(ν)f(ν0),f(ν)=(xcoth(x2)4),x=kBTcmb·\begin{equation} a^{\rm tSZ}_{\nu} = \frac{f(\nu)}{f(\nu_0)},\ f(\nu) = \left(x\coth\left(\frac{x}{2}\right)-4 \right),\ x=\frac{h\nu}{k_\mathrm{B} T_\mathrm{cmb}}\cdot \end{equation}(30)

Thermal SZ × CIB correlation.

Following Like13 the cross-correlation between the thermal SZ and the CIB, tSZ × CIB, is parameterized by a single correlation parameter, ξ, and a fixed template from Addison et al. (2012), (Cν×νtSZ×CIB)=ξ𝒟143tSZ𝒟217CIB×(aνtSZaνCIB+aνtSZaνCIB)×CtSZ×CIB,\begin{eqnarray} \begin{split} \left( C^{{\rm tSZ} \times {\rm CIB}}_{\nu\times\nu'}\right)_{\ell} &= \xi\, \sqrt{{\cal D}^{\rm tSZ}_{143}\, {\cal D}^{\rm CIB}_{217}} \\ &\quad \times \left( a^{\rm tSZ}_{\nu}\, a^{\rm CIB}_{\nu'} + a^{\rm tSZ}_{\nu'}\, a^{\rm CIB}_{\nu} \right) \\ &\quad\times C_{\ell}^{{\rm tSZ} \times {\rm CIB}}, \end{split} \label{eq:tSZCIB} \end{eqnarray}(31)where aνtSZ\hbox{$a^{\rm tSZ}_{\nu}$} is the thermal Sunyaev-Zeldovich spectrum, corrected for the Planck bandpass colour corrections and aνCIB\hbox{$a^{\rm CIB}_{\nu}$} is the CIB spectrum, rescaled at ν = 217 GHz as in the previous paragraphs.

SZ prior.

The kinetic SZ, the thermal SZ, and its correlation with the CIB are not constrained accurately by the Planck data alone. Besides, the tSZ×CIB level is highly correlated with the amplitude of the tSZ. In 2013, we reduced the degeneracy between those parameters and improved their determination by adding the ACT and SPT data. In 2015, we instead impose a Gaussian prior on the tSZ and kSZ amplitudes, inspired by the constraints set by these experiments. From a joint analysis of the Planck 2013 data with those from ACT and SPT, we obtain 𝒟kSZ+1.6𝒟tSZ=(9.5±3)μK2,\begin{equation} {\cal D}^{\rm kSZ} + 1.6 {\cal D}^{\rm tSZ} = (9.5 \pm 3)\,\mu\textrm{K}^2, \end{equation}(32)in excellent agreement with the estimates from Reichardt et al. (2012), once they are rescaled to the Planck frequencies (see Planck Collaboration XIII 2016, for a detailed discussion).

As can be seen in Fig. 17, the kSZ, tSZ, and tSZ×CIB correlations are always dominated by the dust, CIB, and point-source contributions.

3.4. Instrumental modelling

The following sections describe the instrument modelling elements of the model vector, addressing the issues of calibration and beam uncertainties in Sects. 3.4.13.4.3, and describing the noise properties in Sect. 3.4.4. For convenience, Table 10 defines the symbol used for the calibration parameters and the priors later used for exploring them.

3.4.1. Power spectra calibration uncertainties

As in 2013, we allow for a small recalibration of the different frequency power spectra, in order to account for residual uncertainties in the map calibration process. The mixing matrix in the model vector from Eq. (14) can be rewritten as (MZW,ν×νXY)(θinst)=Gν×νXY(θcalib)(MZW,ν×νXY,other)(θother),Gν×νXY(θcalib)=\begin{eqnarray} \left ( {M}^{XY}_{ZW,\nu \times \nu'}\right )_\ell(\theta_{\rm inst}) &=& G^{XY}_{\nu\times\nu'}(\theta_{\rm calib}) \left ( {M}^{XY,\mathrm{other}}_{ZW,\nu \times \nu'}\right )_\ell(\theta_{\rm other}) , \nonumber \\ G^{XY}_{\nu\times\nu'}(\theta_{\rm calib}) &=& \frac{1}{\calibM_{\rm P}^2} \left( \frac{1}{2\sqrt{\calibC^{XX}_{\nu}\calibC^{YY}_{\nu'}}} + \frac{1}{2\sqrt{\calibC^{XX}_{\nu'}\calibC^{YY}_{\nu}}} \right) , \label{eq:caldef} \end{eqnarray}(33)where cνXX\hbox{$\calibC^{XX}_{\nu}$} is the calibration parameter for the XX power spectrum at frequency ν, X being either T or E, and yP is the overall Planck calibration. We ignore the -dependency of the weighting function between the TE and ET spectra at different frequencies that are added to form an effective cross-frequency TE cross-spectrum. As in 2013, we use the TT at 143 GHz as our inter-calibration reference, so that c143TT=1\hbox{$\calibC^{TT}_{143}=1$}.

We further allow for an overall Planck calibration uncertainty, whose variation is constrained by a tight Gaussian prior, yP=1±0.0025.\begin{eqnarray} \calibM_{\rm P} = 1 \pm 0.0025. \end{eqnarray}(34)This prior corresponds to the estimated overall uncertainty, which is discussed in depth in Planck Collaboration I (2016).

The calibration parameters can be degenerate with the foreground parameters, in particular the point sources at high (for TT) and the Galaxy for 217 GHz at low . We thus proceed as in 2013, and measure the calibration refinement parameters on the large scales and on small sky fractions near the Galactic poles. We perform the same estimates on a range of Galactic masks (G20, G30, and G41) restricted to different maximum multipoles (up to = 1500). The fits are performed either by minimizing the scatter between the different frequency spectra, or by using the SMICA algorithm (see Planck Collaboration VI 2014, Sect. 7.3) with a freely varying CMB and generic foreground contribution. For the TT spectra, we obtained in both cases very similar recalibration estimates, from which we extracted the conservative Gaussian priors on recalibration factors, c100TT=0.999±0.001,c217TT=0.995±0.002.\begin{eqnarray} \label{eq:relcal} \calibC^{TT}_{100} &=& 0.999 \pm 0.001 , \\ \calibC^{TT}_{217} &=& 0.995 \pm 0.002 . \end{eqnarray}These are compatible with estimates made at the map level, but on the whole sky; see Planck Collaboration VIII (2016).

3.4.2. Polarization efficiency and angular uncertainty

We now turn to the polarization recalibration case. The signal measured by an imperfect PSB is given by d=G(1+γ)[I+ρ(1+η)(Qcos2(φ+ω)+Usin2(φ+ω))]+n,\begin{equation} d = G(1+\gamma)\left[I + \rho(1+\eta)\left(Q \cos 2(\phi+\omega) + U \sin 2(\phi+\omega) \right) \right] + n , \end{equation}(37)where I, Q, and U are the Stokes parameters; n is the instrumental noise; G, ρ, and φ are the nominal photometric calibration factor, polar efficiency, and direction of polarization of the PSB; and γ, η, and ω are the (small) errors made on each of them (see, e.g., Jones et al. 2007). Due to these errors, the measured cross-power spectra of maps a and b are then contaminated by a spurious signal given by ΔCTT=ΔCTE=ΔCEE=(γa+γb+ηa+ηb2ωa22ωb2)CEE% subequation 9572 0 \begin{eqnarray} \Delta C_\ell^{TT} &=& \left( \gamma_a + \gamma_b \right) C_\ell^{TT}, \label{eq:spurious_Cl_systematics_a} \\ \Delta C_\ell^{\textit{TE}} &=& \left(\gamma_a + \gamma_b + \eta_b - 2 \omega_b^2\right) C_\ell^{\textit{TE}}, \label{eq:spurious_Cl_systematics_b} \\ \Delta C_\ell^{EE} &=& \left(\gamma_a + \gamma_b + \eta_a + \eta_b - 2 \omega_a^2 - 2 \omega_b^2 \right) C_\ell^{EE} \nonumber \\ && \quad + 2 \left(\omega_a^2 + \omega_b^2\right) C_\ell^{BB}, \label{eq:spurious_Cl_systematics_c} \end{eqnarray}where γx, ηx, and ωx, for x = a,b, are the effective instrumental errors for each of the two frequency-averaged maps. Pre-flight measurements of the HFI polarization efficiencies, ρ, had uncertainties | ηx | ≈ 0.3%, while the polarization angle of each PSB is known to | ωx | ≈ 1deg (Rosset et al. 2010). Analysis of the 2015 maps shows the relative photometric calibration of each detector at 100 to 217 GHz to be known to about | γx | = 0.16% at worst, with an absolute orbital dipole calibration of about 0.2%, while analysis of the Crab Nebula observations showed the polarization uncertainties to be consistent with the pre-flight measurements (Planck Collaboration VIII 2016).

Assuming CBB\hbox{$C_\ell^{BB}$} to be negligible, and ignoring ω2 ≪ | η | in Eq. (35), the Gaussian priors on γ and η for each frequency-averaged polarized map would have rms of σγ = 2 × 10-3 and ση = 3 × 10-3. Adding those uncertainties in quadrature, the auto-power spectrum recalibration cνEE\hbox{$\calibC^{EE}_{\nu}$} introduced in Eq. (33) would be given, for an equal-weight combination of nd = 8 polarized detectors, by cνEE=1±2σγ2+ση2nd=1±0.0025.\begin{eqnarray} \calibC^{EE}_{\nu} = 1 \pm 2 \sqrt{ \frac{\sigma_{\gamma}^2 + \sigma_{\eta}^2}{n_\mathrm{d}} } = 1 \pm 0.0025 . \label{eq:calib_EE_prior} \end{eqnarray}(39)The most accurate recalibration factors for TE and EE could therefore be somewhat different from TT. We found, though, that setting the EE recalibration parameter to unity or implementing those priors makes no difference with respect to cosmology; i.e., we recover the same cosmological parameters, with the same uncertainties. Thus, for the baseline explorations, we fixed the EE recalibration parameter to unity, cνEE=1,\begin{equation} \calibC^{EE}_\nu = 1 , \end{equation}(40)and the uncertainty on TE comes only from the TT calibration parameter through Eq. (33).

We also explored the case of much looser priors, and found that best-fit calibration parameters deviate very significantly, and reach values of several percent (between 3% and 12% depending on the frequencies and on whether we fit the EE or TE case). This cannot be due to the instrumental uncertainties embodied in the prior. In the absence of an informative prior, this degree of freedom is used to minimize the differences between frequencies that stem from other effects, not included in the baseline modelling.

The next section introduces one such effect, the temperature-to-polarization leakage, which is due to combining detectors with different beams without accounting for it at the map-making stage (see Sect. 3.4.3). But anticipating the results of the analysis described in Appendix C.3.5, we note that when the calibration and leakage parameters are explored simultaneously without priors, they remain in clear tension with the priors (even if the level of recalibration decreases slightly, by typically 2%, showing the partial degeneracy between the two). In other words, when calibration and leakage parameters are both explored with their respective priors, there is evidence of residual unmodelled systematic effects in polarization – to which we will return.

3.4.3. Beam and transfer function uncertainties

The power spectra from map pairs are corrected by the corresponding effective beam window functions before being confronted with the data model. However, these window functions are not perfectly known, and we now discuss various related sources of errors and uncertainties, the impact of which on the reconstructed Cs is shown in Fig. 21.

Sub-pixel effects.

The first source of error, the so-called “sub-pixel” effect, discussed in detail in Like13, is a result of the Planck scanning strategy and map-making procedure. Scanning along rings with very low nutation levels can result 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 centre of the pixel. This effect, similar to the gravitational lensing of the CMB, has a non-diagonal influence 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. We showed in Like13 that including this effect had little impact on the cosmological parameters measured by Planck.

Masking effect.

A second source of error is the variation, from one sky pixel to another, of the effective beam width, which is averaged over all samples falling in that pixel. While all the HEALPix pixels have the same surface area, their shape – and therefore their moment of inertia (which drives the pixel window function) – depends on location, as shown in Fig. 22, and therefore makes the effective beam window function depend on the pixel mask considered. Of course the actual sampling of the pixels by Planck leads to individual moments of inertia slightly different from the intrinsic values shown here, but spot-check comparisons of this semi-analytical approach used by QuickBeam with numerical simulations of the actual scanning by FEBeCoP showed agreement at the 10-3 level for < 2500 on the resulting pixel window functions for sky coverage varying from 40 to 100%.

In the various Galactic masks used here (Figs. 1213) the contribution of the unmasked pixels to the total effective window function departs from the full-sky average (which is not included in the effective beam window functions), and we therefore expect a different effective transfer function for each mask. We ignored this dependence and mitigated its effect by using transfer functions computed with the Galactic mask G60 which retains an effective sky fraction (including the mask apodization) of fsky = 60%, not too different from the sky fractions fsky between 41 and 70% (see Sect. 3.2.2) used for computing the power spectra.

thumbnail Fig. 21

Contribution of various beam-window-function-related errors and uncertainties to the C relative error. In each panel, the grey histogram shows the relative statistical error on the Planck CMB TT binned power spectrum (for a bin width Δ = 30) divided by 10, while the vertical grey dashes delineate the range < 1800 that is most informative for base ΛCDM. Top: estimation of the error made by ignoring the sub-pixel effects for a fiducial C including the CMB and CIB contributions. Middle: error due to the sky mask, for the Galactic masks used in the TT analysis. Bottom: current beam window function error model, shown at 1σ (solid lines) and 10σ (dotted lines).

Figure 21 compares the impact of these two sources of uncertainty on the stated Planck statistical error bars for Δ = 30. It shows that, for < 1800 where most of the information on ΛCDM lies, the error on the TT power spectra introduced by the sub-pixel effect and by the sky-coverage dependence are less than about 0.1%, and well below the statistical error bars of the binned C. In the range 1800 ≤ ≤ 2500, which helps constrain one-parameter extensions to base ΛCDM (such as Neff), the relative error can reach 0.4% (note as a comparison that the high- ACT experiment states a statistical error of about 3% on the bin 2340 ≤ ≤ 2540, Das et al. 2014). The bottom panel shows the Monte Carlo error model of the beam window functions, which provides negligible (-coupled) uncertainties. Even if this model is somewhat optimistic, since it does not include the effect of the ADC non-linearities and the colour-correction effect of beam measurements on planets (Planck Collaboration VII 2016), we note that even expanding them by a factor of 10 keeps them within the statistical uncertainty of the power spectra.

Modelling the uncertainties.

As in the 2013 analysis, the beam uncertainty eigenmodes were determined from 100 (improved) Monte Carlo (MC) simulations of each planet observation used to measure the scanning beams, then processed through the same QuickBeam pipeline as the nominal beam to determine their effective angular transfer function B(). Thanks to the use of Saturn and Jupiter transits instead of the dimmer Mars used in 2013, the resulting uncertainties are now significantly smaller (Planck Collaboration VII 2016).

For each pair of frequency maps (and frequency-averaged beams) used in the present analysis, a singular-value decomposition (SVD) of the correlation matrix of 100 Monte Carlo based B() realizations was performed over the ranges [0,max] with max = (2000,3000,3000) at (100, 143, 217 GHz), and the five leading modes were kept, as well as their covariance matrix (since the error modes do exhibit Gaussian statistics). We therefore have, for each pair of beams, five -dependent templates, each associated with a Gaussian amplitude centred on 0, and a covariance matrix coupling all of them.

thumbnail Fig. 22

Map of the relative variations of the trace of the HEALPix pixel moment of inertia tensor at Nside = 2048 in Galactic coordinates.

Including the beam uncertainties in the mixing matrix of Eq. (14) gives (MZW,ν×νXY)(θinst)=(MZW,ν×νXY,other)(θother)(ΔWν×νZW)(θbeam),(ΔWν×νZW)(θbeam)=expi=152θν×νZW,i(Eν×νZW,i),\begin{eqnarray} \left ( M^{XY}_{ZW,\nu \times \nu'}\right )_\ell(\theta_{\rm inst}) &=& \left ( M^{XY,\mathrm{other}}_{ZW,\nu \times \nu'}\right )_\ell(\theta_{\rm other})\ \left( \Delta W^{ZW}_{\nu \times \nu'} \right )_\ell(\theta_{\rm beam}) , \nonumber \\ \left( \Delta W^{ZW}_{\nu \times \nu'} \right )_\ell(\theta_{\rm beam}) &= &\exp{\sum_{i=1}^5 2\, \theta^{ZW,i}_{\nu\times\nu'} \left (E^{ZW,i}_{\nu\times\nu'}\right )_\ell} , \end{eqnarray}(41)where (ΔWν×νZW)(θbeam)\hbox{$\left( \Delta W^{ZW}_{\nu \times \nu'} \right )_\ell(\theta_{\rm beam})$} stands for the beam error built from the eigenmodes Eν×νZW,i)(\hbox{$\left (E^{ZW,i}_{\nu\times\nu'}\right )_\ell$}. The quadratic sum of the beam eigenmodes is shown in Fig. 21. This is much smaller (less than a percent) than the combined TT spectrum error bars. This contrasts with the 2013 case where the beam uncertainties were greater; for instance, for the 100, 143, and 217 GHz channel maps, the rms of the W() = B()2 uncertainties at = 1000 dropped from (61,23,20) × 10-4 to (2.2,0.84,0.81) × 10-4, respectively. The fact that beam uncertainties are sub-dominant in the total error budget is even more pronounced in polarization, where noise is higher. We use the beam modes computed from temperature data, combined with appropriate weights when used as parameters affecting the TE and EE spectra.

As in 2013, instead of including the beam error in the vector model, we include its contribution to the covariance matrix, linearizing the vector model so that (Cν×νXY)(θ)=(Cν×νXY)(θ,θbeam=0)+(ΔWν×νZW)(θbeam)(Cν×νXY),\begin{equation} \left ( C^{XY}_{\nu \times \nu'} \right )_\ell(\theta) = \left ( C^{XY}_{\nu \times \nu'} \right )_\ell(\theta,\theta_{\rm beam}=0) + \left( \Delta W^{ZW}_{\nu \times \nu'} \right )_\ell(\theta_{\rm beam}) \left (C^{XY}_{\nu \times \nu'} \right )^*_\ell , \end{equation}(42)where Cν×νXY)(\hbox{$\left (C^{XY}_{\nu \times \nu'} \right )^*_\ell$} is the fiducial spectrum XY for the pair of frequencies ν × ν obtained using the best cosmological and foreground model. We can then marginalize over the beam uncertainty, enlarging the covariance matrix to obtain Cbeammarg.=C+CΔWΔWTCT,\begin{equation} \tens{C}_{\mathrm{beam\ marg.}} = \tens{C} + \vec{C}^{*} \left < \Delta \vec{W}\Delta \vec{W}^\tens{T}\right > \vec{C}^{*\tens{T}} , \end{equation}(43)where ΔWΔWT\hbox{$\left < \Delta\vec{W}\Delta\vec{W}^\tens{T}\right >$} is the Monte Carlo based covariance matrix, restricted to its first five eigenmodes.

In 2013, beam errors were marginalized for all the modes except the two greatest of the 100 × 100 spectrum. In the present release we instead marginalize over all modes in TT, TE, and EE. We also performed a test in which we estimated the amplitudes for all of the first five beam eigenmodes in TT, TE, and EE, and found no indication of any beam error contribution (see Sect. 4.1.3 and Fig. 35).

Temperature-to-polarization leakage.

Polarization measurements are differential by nature. Therefore any unaccounted discrepancy in combining polarized detectors can create some leakage from temperature to polarization (Hu et al. 2003). Sources of such discrepancies in the current HFI processing include, but are not limited to: differences in the scanning beams that are ignored during the map-making; differences in the noise level, because of the individual inverse noise weighting used in HFI; and differences in the number of valid samples.

For this release, we did not attempt to model and remove a priori the form and amplitude of this coupling between the measured TT, TE, and EE spectra; we rather estimate the residual effect by fitting a posteriori in the likelihood some flexible template of this coupling, parameterized by some new nuisance parameters that we now describe.

The temperature-to-polarization leakage due to beam mismatch is assumed to affect the spherical harmonic coefficients via aℓmTaℓmT,aℓmE% subequation 10071 0 \begin{eqnarray} a_{\ell m}^T &\longrightarrow& a_{\ell m}^T , \\[3mm] \label{eq:bleak_alm_t} a_{\ell m}^E &\longrightarrow& a_{\ell m}^E + \varepsilon({\ell}) a_{\ell m}^T , \label{eq:bleak_alm_e} \end{eqnarray}and, for each map, the spurious polarization power spectrum CXYmaℓmXaℓmY/(2+1)\hbox{$C_\ell^{XY} \equiv \sum_m a^X_{\ell m} a^{Y*}_{\ell m} / (2\ell+1)$} is modelled as ΔCTE=ΔCEE=% subequation 10086 0 \begin{eqnarray} \Delta C_{\ell}^{\textit{TE}} &=& \varepsilon({\ell})\, C_{\ell}^{TT}, \label{eq:bleak_cl_a} \\[3mm] \Delta C_{\ell}^{EE} &=& \varepsilon^2({\ell})\, C_{\ell}^{TT} + 2 \varepsilon({\ell})\, C_{\ell}^{\textit{TE}}. \label{eq:bleak_cl_b} \end{eqnarray}Here ε is a polynomial in multipole determined by the effective beam of the detector-assembly measuring the polarized signal. Considering an effective beam map b(nˆ)\hbox{$b(\hat{\vec n})$} (rotated so that it is centred on the north pole), its spherical harmonic coefficients are defined as bℓmdnˆb(nˆ)Yℓm(nˆ)\hbox{$b_{\ell m} \equiv \int {\rm d}{\hat{\vec n}}\, b(\hat{\vec n})\, Y^*_{\ell m}(\hat{\vec n})$}. As a consequence of the Planck scanning strategy, pixels are visited approximately every six months, with a rotation of the focal plane by 180deg, and we expect bℓm to be dominated by even values of m, and especially the modes m = 2 and 4, which describe the beam ellipticity. As noted by, e.g., Souradeep & Ratra (2001) for elliptical Gaussian beams, the Planck-HFI beams for a detector d obey bℓm(d)βm(d)mb0(d).\begin{equation} b_{\ell m}^{(d)} \simeq \beta_m^{(d)} \ell^m b_{\ell 0}^{(d)} . \label{eq:bleak_lpower} \end{equation}(46)We therefore fit the spectra using a fourth-order polynomial ε()=ε0+ε22+ε44,\begin{equation} \varepsilon(\ell) = \varepsilon_0 + \varepsilon_2 \ell^2 + \varepsilon_4 \ell^4, \label{eq:bleak_pol} \end{equation}(47)treating the coefficients ε0, ε2, and ε4 as nuisance parameters in the MCMC analysis. Tests performed on detailed simulations of Planck observations with known mismatched beams have shown that Eqs. (42) and (47) describe the power leakage due to beam mismatch with an accuracy of about 20% in the range 1002000.

thumbnail Fig. 23

Best fit of the power spectrum leakage due to the beam mismatch for TE (Eq. (45a), upper panel) and EE (Eq. (45b), lower panel). In each case, we show the correction for individual cross-spectra (coloured thin lines) and the co-added correction (black line). The individual cross-spectra corrections are only shown in the range of multipoles where the data from each particular pair is used. The individual correction can be much higher than the co-added correction. The co-added correction is dominated by the best S/N pair for each multipole. For example, up to = 500, the TE co-added correction is dominated by the 100 × 143 contribution. The grey dashed lines show the TE and EE best-fit spectra rescaled by a factor of 20, to give an idea of the location of the model peaks.

The equations above suggest that the same polynomial ε can describe the contamination of the TE and EE spectra for a given pair of detector sets. But in the current Plik analysis, the TE cross-spectrum of two different maps a and b is the inverse-variance-weighted average of the cross-spectra TaEb and TbEa, while EE is simply EaEb. In addition, the temperature maps include the signal from SWBs, which is obviously not the case for the E maps. We therefore describe the TE and EE corrections by different ε parameters. Similarly, we treated the parameters for the EE cross-frequency spectra as being uncorrelated with the parameters for the auto-frequency ones.

The leakage is driven by the discrepancy between the individual effective beams bℓm(d)\hbox{$b_{\ell m}^{(d)}$} making up a detector assembly, coupled with the details of the scanning strategy and relative weight of each detector. If we assumed a perfect knowledge of the beams, precise – but not necessarily accurate – numerical predictions of the leakage would be possible. However, we preferred to adopt a more conservative approach in which the leakage was free to vary over a range wide enough to enclose the true value. On the other hand, in order to limit the unphysical range of variations permitted by so many nuisance parameters, we need priors on the εm terms used in the Monte Carlo explorations. We assume Gaussian distributions of zero mean with a standard deviation σm representative of the dispersion found in simulations of the effect with realistic instrumental parameters. We found σ0 = 1 × 10-5, σ2 = 1.25 × 10-8, and σ4 = 2.7 × 10-15. This procedure ignores correlations between terms of different m, and is therefore likely substantially too permissive.

Another way of deriving the beam leakage would be to use a cosmological prior, i.e., by finding the best fit when holding the cosmological parameters fixed at their best-fit values for base ΛCDM. Figure 23 shows the result of this procedure for the cross-frequency pairs. The figure also shows the implied correction for the co-added spectra. This correction is dominated by the pair with the highest S/N at each multipole. The fact that different sets are used in different -ranges leads to discontinuities in the correction template of the co-added spectrum. As can be seen in the figure, the co-added beam-leakage correction, of order μK2, is much smaller than the individual corrections, which partially compensate each other on average (but improve the agreement between the individual polarized cross-frequency spectra).

It is shown in Appendix C.3.5 that neither procedure is fully satisfactory. The cosmological prior leads to nuisance parameters that vastly exceed the values allowed by the physical priors, and the physical priors are clearly overly permissive (leaving the cosmological parameters unchanged but with doubled error bars for some parameters). In any case, the agreement between the different cross-spectra remains much poorer in polarization than in temperature (see Sect. 4.4, Fig. 40, and Appendix C.3.5); they present oscillatory features similar to the ones produced by our beam leakage model, but the model is clearly not sufficient. For lack of a completely satisfactory global instrumental model, this correction is only illustrative and it is not used in the baseline likelihood.

3.4.4. Noise modelling

To predict the variance of the empirical power spectra, we need to model the noise properties of all maps used in the construction of the likelihood. As described in detail in Planck Collaboration VII (2016) and Planck Collaboration VIII (2016), the Planck HFI maps have complicated noise properties, with noise levels varying spatially and with correlations between neighbouring pixels along the scanning direction.

thumbnail Fig. 24

Deviations from a white noise power spectrum induced by noise correlations. We show half-ring difference power spectra for 100 GHz half-mission 1 maps (blue lines) of Stokes parameters I (top panel), Q (middle panel), and U (bottom panel). The best-fitting analytical model of the form Eq. (48) is over-plotted in red.

thumbnail Fig. 25

Difference between auto and cross-spectra for the 100 GHz half-mission maps, divided by the noise estimate from half-ring difference maps (blue and green lines). Noise estimates derived from half-ring difference maps are biased low. We fit the average of both half-mission curves (black line) with a power law model (red line). The analysis procedure is applied to the Stokes parameter maps I, Q, and U (top to bottom). All data power spectra are smoothed.

For each channel, full-resolution noise variance maps are constructed during the map-making process (Planck Collaboration VIII 2016). They provide an approximation to the diagonal elements of the true npix × npix noise covariance matrix for Stokes parameters I (temperature only), or I, Q, and U (temperature and polarization). While it is possible to capture the anisotropic nature of the noise variance with these objects, noise correlations between pixels remain unmodelled. To include deviations from a white-noise power spectrum, we therefore make use of half-ring difference maps. Choosing the 100 GHz map of the first half-mission as an example, we show the scalar (spin-0) power spectra of the three temperature and polarization maps in Fig. 24, rescaled by arbitrary constants. We find that the logarithm of the HFI noise power spectra as given by the half-ring difference maps can be accurately parameterized using a fourth-order polynomial with an additional logarithmic term, log(CHRD)=i=04αii+α5log(+α6).\begin{equation} \log(C^{\mathrm{HRD}}_{\ell}) = \sum_{i = 0}^{4} \alpha_i \, \ell^{i} + \alpha_5 \log(\ell + \alpha_6) . \label{eq:plik_hrd_noise_fit} \end{equation}(48)Since low-frequency noise and processing steps like deglitching leave residual correlations between both half-ring maps, noise estimates derived from their difference are biased low, at the percent level at high- (where it was first detected and understood, see Planck Collaboration VI 2014). We correct for this effect by comparing the difference of auto-power-spectra and cross-spectra (assumed to be free of noise bias) at a given frequency with the noise estimates obtained from half-ring difference maps. As shown in Fig. 25, we use a a power-law model with free spectral index to fit the average of the ratios of the first and second half-mission results to the half-ring difference spectrum, using the average to nullify chance correlations between signal and noise: Cbias=α0α1+α2.\begin{equation} C^{\mathrm{bias}}_{\ell} = \alpha_0 \, \ell^{\alpha_1} + \alpha_2 . \label{eq:plik_hrd_noise_bias} \end{equation}(49)At a multipole moment of = 1000, we obtain correction factors for the temperature noise estimate obtained from half-ring difference maps of 9%, 10%, and 9% at 100, 143, and 217 GHz, respectively.

In summary, our HFI noise model is obtained as follows. For each map, we capture the anisotropic nature of the noise amplitude by using the diagonal elements of the pixel-space noise covariance matrix. The corresponding white-noise power spectrum is then modulated in harmonic space using the product of the two smooth fitting functions given in Eqs. (48) and (49).

thumbnail Fig. 26

Correlated noise model. In grey are shown the cross-detector TT spectra of the half-ring difference maps. The black line show the same, smoothed by a Δ = 200 sliding average, while the blue data points are a Δ = 100 binned version of the grey line. Error bars simply reflect the scatter in each bin. The green line is the spline-smoothed version of the data that we use as our correlated noise template.

Correlated noise between detectors.

If there is some correlation between the noise in the different cuts in our data, the trick of only forming effective frequency-pair power spectra from cross-spectra to avoid the noise biases fails. In 2013, we evaluated the amplitude of such correlated noise between different detsets. The correlation, if any, was found to be small, and we estimated its effect on the cosmological parameter fits to be negligible. As stated in Sect. 3.2.1, the situation is different for the 2015 data. Indeed, we now detect a small but significant correlated noise contribution between the detsets. This is the reason we change our choice of data to estimate the cross-spectra, from detsets to half-mission maps. The correlated noise appears to be much less significant in the latter.

To estimate the amount of correlated noise in the data, we measured the cross-spectra between the half-ring difference maps of all the individual detsets. The cross-spectra are then summed using the same inverse-variance weighting that we used in 2013 to form the effective frequency-pair spectra. Figure 26 shows the spectra for each frequency pair. All of these deviate significantly from zero. We build an effective correlated noise template by fitting a smoothing spline on a Δ = 200 sliding average of the data. Given the noise level in polarization, we did not investigate the possible contribution of correlated noise in EE and TE.

Section 4.1.1 shows that when these correlated noise templates are used, the results of the detsets likelihood are in excellent agreement with those based on the baseline, half-mission one.

3.5. Covariance matrix structure

The construction of a Gaussian approximation to the likelihood function requires building covariance matrices for the pseudo-power spectra. Mathematically exact expressions exist, but they are prohibitively expensive to calculate numerically at Planck resolution (Wandelt et al. 2001); we thus follow the approach taken in Like13 and make use of analytical approximations (Hansen et al. 2002; Hinshaw et al. 2003; Efstathiou 2004; Challinor & Chon 2005).

For our baseline likelihood, we calculate covariance matrices for all 45 unique detector combinations that can be formed out of the six frequency-averaged half-mission maps at 100, 143, and 217 GHz. To do so, we assume a fiducial power spectrum that includes the data variance induced by the CMB and all foreground components described in Sect. 3.3; this variance is computed assuming these components are Gaussian-distributed. The effect of this approximation regarding Galactic foregrounds is tested by means of simulations in Sect. 3.6. The fiducial model is taken from the best-fit cosmological and foreground parameters; since they only become available after a full exploration of the likelihood, we iteratively refine our initial guess. As discussed in Sect. 3.1, the data vector used in the likelihood function of Eq. (13) is constructed from frequency-averaged power spectra. Following Like13, for each polarization combination, we therefore build averaged covariance matrices for the four frequencies ν1,ν2,ν3,ν4, Var(XYν1,ν2,ZWν3,ν4)=(i,j)(ν1,ν2)(p,q)(ν3,ν4)wXYi,jwZWp,q\begin{eqnarray} \mathrm{Var}(\hat{C}_{\ell}^{XY \ \nu_1, \nu_2}, \hat{C}_{\ellp}^{ZW \ \nu_3, \nu_4}) &=& \sum_{\substack{(i, j) \in (\nu_1, \nu_2) \\ (p, q) \in (\nu_3, \nu_4)}} w^{XY \ i,j}_{\ell} w^{ZW \ p, q}_{\ellp} \nonumber\\ &&\quad \times \mathrm{Var}\left(\hat{C}_{\ell}^{XY \ i, j}, \hat{C}_{\ellp}^{ZW \ p, q}\right) , \label{eq:hil_cl_cov_avg} \end{eqnarray}(50)where X,Y,Z,W ∈ { T,E }, and wXYi,j is the inverse-variance weight for the combination (i,j), computed from wXYi,j1/Var(XYi,j,XYi,j),\begin{equation} w^{XY \ i,j}_{\ell} \propto 1/\mathrm{Var}\left(\hat{C}_{\ell}^{XY \ i, j}, \hat{C}_{\ell}^{XY \ i, j}\right) , \end{equation}(51)and normalized to unity. For the averaged XY = TE covariance (and likewise for ZW = TE), the sum in Eq. (50) must be taken over the additional permutation XY = ET. That is, the two cases where the temperature map of channel i is correlated with the polarization map of channel j and vice versa are combined into a single frequency-averaged covariance matrix. These matrices are then combined to form the full covariance used in the likelihood, C=(CTTTTCTTEECTTTECEETTCEEEECEETECTETTCTEEECTETE),\begin{equation} \tens{C} = \begin{pmatrix} C^{TTTT} & C^{TTEE} & C^{TTTE} \\ C^{EETT} & C^{EEEE} & C^{EETE} \\ C^{TETT} & C^{TEEE} & C^{TETE} \end{pmatrix} , \end{equation}(52)where the individual polarization blocks are constructed from the frequency-averaged covariance matrices of Eq. (50) (Like13).

Table 13

Shifts of parameters over 300 TT simulations.

Appendix C.1.1 provides a summary of the equations used to compute temperature and polarization covariance matrices and presents a validation of the implementation through direct simulations. Let us note that, for the approximations used in the analytical computation of the covariance matrix to be precise, the mask power spectra have to decrease quickly with multipole moment ; this requirement gives rise to the apodization scheme discussed in Sect. 3.2.2. In the presence of a point-source mask, however, the condition may no longer be fulfilled, reducing the accuracy of the approximations assumed in the calculation of the covariance matrices. We discuss in Appendix C.1.4 the heuristic correction we developed to restore the accuracy, which is based on direct simulations of the effect.

3.6. FFP8 simulations

In order to validate the overall implementation and our approximations, we generated 300 simulated HFI half-mission map sets in the frequency range 100 to 217 GHz, which we analysed like the real data. For the CMB, we created realizations of the ΛCDM model with the best-fit parameters obtained in this paper. After convolving the CMB maps with beam and pixel window functions, we superimposed CIB, dust, and noise realizations from the FFP8 simulations (Planck Collaboration XII 2016) that capture both the correlation structure and anisotropy of foregrounds and noise. We then computed power spectra using the set of frequency-dependent masks described in Sect. 3.2.2 and created the corresponding Plik TT likelihood. We modified the shape of the foreground spectra to fit the FFP8 simulations, but kept the parameterization used on the data. In the case of dust, we used priors similar to those used on data. Furthermore, in the following the dust amplitude parameter is named gal545ν×ν\hbox{$\mathrm{gal}^{\nu\times\nu'}_{545}$}. We then ran an MCMC sampler to derive the cosmological and foreground parameters posterior distributions for all dataset realizations.

thumbnail Fig. 27

Plik parameter results on 300 simulations for the six baseline cosmological parameters, as well as the FFP8 CIB and Galactic dust amplitudes. The simulations include quite realistic CMB, noise, and foregrounds (see text). The distributions of inferred posterior mean parameters are centred around their input values with the expected scatter. Indeed the dotted red lines show the best-fit Gaussian for each distribution, with a mean shift, Δμ, and a departure Δσ from unit standard deviation given in the legend; both are close to zero. These best fits are thus very close to Gaussian distributions with zero shift and unit variance, which are displayed for reference as black lines. The legend gives the numerical value of Δμ and Δσ, as well as the p-values of a Kolmogorov-Smirnov test of the histograms against a Gaussian distribution shifted from zero by Δμ and with standard deviation shifted from unity by Δσ. This confirms that the distributions are consistent with Gaussian distributions with zero mean and unit standard deviation, with a small offset of the mean.

For each simulation, we computed the shift of the derived posterior mean parameters with respect to the input cosmology, normalized by their posterior widths σpost. When a Gaussian prior with standard deviation σprior is used, we rescale σpost by [1σpost2/σprior2]1/2\hbox{$[1-\sigma_{\mathrm{post}}^2 / \sigma_{\mathrm{prior}}^2]^{1/2}$}; this is the case for τ and for the Galactic dust amplitudes gal545ν\hbox{$\mathrm{gal}^\nu_{545}$} in the four cross-frequency channels used. In Fig. 27, we show histograms of the shifts we found for all 300 simulations for the six baseline cosmological parameters, as well as the FFP8 CIB and galactic dust amplitudes. As shown in the figure, we recover the input parameters with little bias and a scatter of the normalized parameter shifts around unity. The p-values of the Kolmogorov-Smirnov test that we ran are given in the legend and we do not detect significant departures from normality. The average reduced χ2 for the histograms of Fig. 27 is equal to 1.02.

Table 13 (second column) compiles the average shifts of Fig. 27, but in order to gauge whether they are as small as expected for this number of simulations (assuming no bias), the shifts are expressed in units of the posterior width rescaled by 1/300\hbox{$1{/}\!\sqrt{300}$}. We note that the shift of the average is above one (scaled)σ in three cases out of a total of 11 parameters (68% of the Δs would be expected to lie within 1σ if the parameters were uncorrelated), with θ, ns, and gal545217\hbox{$\mathrm{gal}_{545}^{217}$} at the 1.7, 2.0, and 1.5 (scaled)σ level, respectively.

Before proceeding, let us note that an estimate (third column) of these shifts is obtained by simply computing the shift from a single likelihood using as input the average spectra of the 300 simulations. This effectively reduces cosmic variance and noise amplitude by a factor 300\hbox{$\sqrt{300}$} and, more importantly, it decreases the cost and length of the overall computation, enabling additional tests. These shift estimates are noted rA. The table shows that significant improvement in the determination of ns is obtained by removing low- multipoles. Indeed, Cols. 4 and 5 of Table 13 show the variation of the shift when the min of the high- likelihood is increased from 30 to 65 and 100. The shift in ns is decreased by a factor two, while the decrease in the number of bins per cross-frequency spectrum is only reduced from 199 to 185 (having little impact on the size of the covariance matrix of cosmological parameters).

These changes with min therefore trace the small biases back to the lowest- bins. It suggests that the Gaussian approximation used in the high- likelihood starts to become mildly inaccurate at = 30. Indeed, even if noticeable, this effect would contribute at most a 0.11σ bias on ns. This is further confirmed by the lack of a detectable effect found in Sect. 5.1 when varying the hybridization scale in TT between Commander and Plik . However, the exclusion of low- information degrades our ability to accurately reconstruct the foreground amplitudes ACIB217\hbox{$A_{\mathrm{CIB}}^{217}$}, gal545143217\hbox{$\mathrm{gal}_{545}^{143-217}$}, and gal545217\hbox{$\mathrm{gal}_{545}^{217}$}. Indeed, the dust spectral amplitudes in the 143 × 217 and 217 × 217 channels are highest at low multipoles, and the CIB spectrum in the range 30 ≤ ≤ 100 also adds substantial information.

In spite of this low- trade-off between an accurate determination of ns on the one hand and ACIB217\hbox{$A_{\mathrm{CIB}}^{217}$}, gal545143217\hbox{$\mathrm{gal}_{545}^{143-217}$}, and gal545217\hbox{$\mathrm{gal}_{545}^{217}$} on the other, we can conclude that the Plik implementation is behaving as expected and can be used for actual data analysis.

Appendix C.2 extends this conclusion to the joint Plik TT, EE, TE likelihood case.

3.7. End-to-end simulations

Table 14

End-to-end parameter shifts for a single realization of CMB and foregrounds, along with five different noise realizations. Shifts are computed with respect to those obtained without noise and with instrumental effects turned off.

Table 15

End-to-end parameter shifts for four different CMB realizations but comprising four pairs of realizations with the same noise realization with respect to those obtained without noise and with instrumental effects turned off.

While the previous section validated our methodology, our approximations, and the overall implementation, this does not yet give the sensitivity to residual systematic uncertainties undetected by data consistency checks. These are by their very nature very much more difficult to address realistically, since, when an effect is detected and sufficiently well understood, it can be modelled and is corrected for, in general at the TOI-processing stage; only the uncertainty of the correction needs to be addressed. Still, HFI has developed a complete model of the instrument which contains all identified systematic effects and enables realistic simulation of the instrumental response. We have therefore generated a number of full-mission time streams which we have then processed with the DPC TOI processing pipeline in order to create map datasets as close to instrumental reality as we can in order to assess the possible impact of low-level residual instrumental systematics, the effects of which might have remained undetected otherwise.

In this section, we report on the shifts in the values of the cosmological and foreground parameters induced by these specific residual systematic effects, comparing the results of a TT likelihood analysis for two overlapping sets of five simulations:

  • 1.

    five simulations of maps at 100 GHz,143 GHz and 217 GHz, for asingle realization of the CMB and of the foregrounds but for fivedifferent realizations of the noise,

  • 2.

    five simulations of maps at 100 GHz, 143 GHz and 217 GHz, composed of four CMB realizations, two noise realizations, a single realization of the foregrounds, but forming four pairs of realizations having the same noise but different CMB.

These simulations sum to a total of eight distinct simulations and are numbered from 1 to 8 in Tables 14 and 15. To be more explicit, in the second set, among simulations numbered 4 to 8, simulations 4, 7 and 8 have different CMB, but the same noise as each other. Simulations 5 and 6 have different CMB, and the same noise as each other, but different from simulations 4, 7 and 8. Realizations 4 and 5, having the same CMB but different noise, are common to the two sets of five realizations.

Each of these have been performed twice, with the end-to-end (instrument plus TOI processing) pipeline and noise contribution switched either on or off. End-to-end simulations are computationally very costly, (typically a week for each simulated mission dataset) and hence only a few realizations were generated).

As explained in Sect. 5.4 of Planck Collaboration VII (2016), the end-to-end simulations are created by feeding the TOI processing pipeline with simulated data to evaluate and characterize the overall transfer function and the respective contribution of each individual effect on the determination of the cosmological parameters. Simulated TOIs are produced by applying the real mission scanning strategy to a realistic input sky specified by the Planck Sky Model (PSM; Delabrouille et al. 2013) containing a lensed CMB realization, galactic diffuse foregrounds, and the dipole components. To this sky-scanned TOI, we add a white-noise component, representing the phonon and photon noises. The very low-temporal-frequency thermal drift seen in the real data is also added to the TOI. The noisy sky TOI is then convolved with the appropriate bolometer transfer functions. Another white-noise component, representing Johnson noise and read-out noise, is also added. Simulated cosmic rays using the measured glitch rates, amplitudes, and shapes are added to the TOI. This TOI is interpolated to the electronic HFI fast-sampling frequency. It is then converted from analogue to digital using a simulated non-linear analogue-to-digital converter (ADC). Identified 4 K cooler spectral lines are added to the TOI. Both effects (ADC and 4 K lines) are derived from the measured in-flight behaviour. The TOI finally goes through the data compression/decompression algorithm used for communication between the Planck satellite and Earth. The simulated TOI is then processed in the same way as the real mission data for cleaning and systematic error removal, calibration, destriping, and map-making.

Some limitations of the current end-to-end approach follow. No pointing error is included, although previous (dedicated) simulations suggest that this has negligible effect. In addition, this effect was included in the dedicated simulations performed to assess the precision of the beam recovery procedure. The first step of the TOI processing is to correct the ADC non-linearity (ADC NL). For the flight data, the ADC NL was determined by using HFI’s measured signal at the end of the HFI mission, with the instrument’s cooling system switched off and an instrument temperature equal to 4 K. This determination relied on supposing the signal to be perfect white noise and therefore to correspond to the distortions brought in by the ADC. In the current implementation, we assume perfect knowledge of ADC NL and 4 K lines. This is of course not true for the real data and future end-to-simulations, accompanying Planck’s next data release, will improve our model of this effect. After ADC NL correction, the signal is converted to volts. Deglitching is then performed by flagging glitch heads and by using glitch tail time-lines. This enables the creation of the thermal baseline which is used for signal demodulation. The thermal baseline and glitch tails are subtracted, the signal is converted to watts, and the 4 K lines are removed. The resulting signal is then deconvolved by the bolometer transfer functions. We do not include uncertainties in the glitch tail shape used in the deglitching procedure, i.e., the templates are the same for the simulations and the processing; but here again, previous studies suggest any difference is a small effect.

The analysis of these sets of end-to-end simulations, and of their counterparts for which all instrumental effects are turned off, is performed similarly to that of the simulations described in Sect. 3.6. Angular power spectra for all cross-half-missions and for all frequency combinations are computed using the Planck masks described in Sect. 3.2.2 and with the appropriate beam functions. Noise levels are evaluated as described in Sect. 3.4.4. Templates for galactic foregrounds (CO, free-free, synchrotron, thermal and spinning dust), the kinetic and thermal SZ effects, the cosmic infrared background, and radio and IR point sources are constructed based on the PSM input foreground maps. The covariance matrix is computed with the method outlined in Appendix C.1 with the aforementioned input CMB power spectrum, input foreground spectra, noise levels, beam functions and masks.

thumbnail Fig. 28

Plik 2015 co-added TT, TE, and EE spectra. The blue points are for bins of Δ = 30, while the grey points are unbinned. The lower panels show the residuals with respect to the best fit Plik TT+tauprior ΛCDM model. The yellow lines show the 68% unbinned error bars. For TE and EE, we also show the best-fit beam-leakage correction (green line; see text and Fig. 23).

thumbnail Fig. 29

Zoom in to various ranges of the HM co-added power spectra, together with the Plik TT+tauprior ΛCDM best-fit model (red line). We show the TT (top), TE (centre) and EE (bottom) power spectra. The lower panels in each plot show the residuals with respect to the best-fit model.

All sets of power spectra and the inverse covariance matrix are then binned and used in the likelihood analysis performed using an MCMC sampler together with Plik and PICO in order to determine the best fit cosmological parameters. The shifts in cosmological parameter values induced by the imperfect correction of instrumental effets by the TOI processing pipeline are then computed for the end-to-end simulations with respect to those obtained for the simulations without noise and with instrumental effects turned off, normalized by the end-to-end simulations’ posterior widths. Comparing shifts computed in this way cancels out cosmic variance and chance correlations between the CMB and the foregrounds and are thus fully attributable to the instrument and to the noise, which cannot be disentangled, as well as to CMB-noise chance correlations. That is, those shifts probe directly the scatter and possible biases induced by residual systematics effects.

The mean and median shifts for the five simulations with a single CMB realization and a single foreground realization, but different noise realizations, are given in Table 14. In order to verify that these shifts are within expectations, we computed the shifts in cosmological parameters for 100 FFP8 simulations, each with identical CMB signal but different FFP8 noise, with respect to the cosmological parameters obtained for the CMB only, normalized by their posterior widths. The standard deviations of the resulting distributions are given in the column labelled “σFFP8” of Table 14 and can be compared with the shifts obtained for the five end-to-end simulations. All shifts are within 1σ of the shifts expected from FFP8. In addition, there is no indication of any detectable bias. All shifts are thus compatible with scatter introduced by noise.

The shifts for the five realizations with four different CMB realizations, the same foregrounds, but comprising four pairs with the same noise realization, are given in Cols. 4 to 8 of Table 15. As mentioned at the beginning of this section, realizations numbered 4 and 5 are common to the sets of Tables 14 and 15. In the columns labelled “Δ5−6” to “Δ7−8”, we computed the absolute differences in the shifts within pairs of realizations having different CMB but the same foreground and noise realizations. We compare these differences to the standard deviations of the distributions of cosmological parameter shifts of 100 FFP8 simulations, varying the CMB but keeping the same FFP8 noise realization, with respect to the cosmological parameters of the corresponding CMB but without noise (column labelled “σFFP8”). These distributions quantify the impact of CMB-noise correlations on the determination of the cosmological parameters. The Table shows that among all Δ’s, 11 are within 1σFFP8, 7 are within 1 to 2σFFP8, 3 are within 2 to 3σFFP8, 2 are within 3 to 4σFFP8. 50% of the differences are within 1σ and 78% within 2σ. At the very worst, taking the example of Δ7−8 a cosmological parameter (θ) moves a total of 4σ, from −3σ to 1σ in units of σFFP8 when the CMB is changed but the noise is left the same. This is rare but can be expected in a few percent of simulations. As in the case of the shifts listed in Table 14, there is thus no detectable bias, with all shifts compatible with those expected from FFP8.

In summary, we have detected no sign as yet of systematic biases of the cosmological parameters due to known low-level instrumental effects as corrected by the current HFI TOI processing pipeline. An increase in the significance of these tests is left for further work once the simulation chain is further optimized for more massive numerical work.

thumbnail Fig. 30

Residuals of the co-added CMB TT power spectra, with respect to the Plik TT+tauprior best-fit model, in units of standard deviation. The three coloured bands (from the centre, yellow, orange, and red) represent the ± 1, ± 2, and ± 3σ regions.

thumbnail Fig. 31

Inter-frequency foreground-cleaned TT power spectra differences, in μK2. Each of the sub-panels shows the difference, after foreground subtraction, between pairs of frequency power spectra (the spectrum named on the vertical axis minus the one named on the horizontal axis), in units of standard deviation. The coloured bands identify deviations that are smaller than one (yellow), two (orange), or three (red) standard deviations. We show the differences for both the HM power spectra (blue points) and the DS power spectra (light blue points) after correlated noise correction. Figure 41 displays the same quantities for the TE and EE spectra.

3.8. High-multipole reference results

This section describes the results obtained using the baseline Plik likelihood, in combination with a prior on the optical depth to reionization, τ = 0.07 ± 0.02 (referred to, in TT, as Plik TT+tauprior). The robustness and validation of these results (presented in Sect. 4) can therefore be assessed independently of any potential low- anomaly, or hybridization issues. The full low- + high- likelihood will be discussed in Sect. 5.

Figure 28 shows the high- co-added CMB spectra in TT, TE, and EE, and their residuals with respect to the best-fit ΛCDM model in TT (red line), both -by- (grey points) and binned (blue circles). The blue error bars per bin are derived from the diagonal of the covariance matrix computed with the best-fit CMB as fiducial model. The bottom sub-panels with residuals also show (yellow lines) the diagonal of the -by- covariance matrix, which may be compared to the dispersion of the individual determinations. Parenthetically, it provides graphical evidence that TT is dominated by cosmic variance through ≈ 1600, while TE is cosmic-variance dominated at ≲ 160 and ≈ 260−460. The jumps in the polarization diagonal-covariance error-bars come from the variable ranges retained at different frequencies, which therefore vary the amount of data included discontinuously with . Figure 29 zooms in to five adjacent ranges on the co-added spectra to allow close inspection of the data distribution around the model.

More quantitatively, Table 16 shows the χ2 values with respect to the ΛCDM best fit to the Plik TT+tauprior data combination for the unbinned CMB co-added power spectra (obtained as described in Appendix C.4). The TT spectrum has a reduced χ2 of 1.03 for 2479 degrees of freedom, corresponding to a probability to exceed (PTE) of 17.2%; the base ΛCDM model is therefore in agreement with the co-added data. The best-fit ΛCDM model in TT also provides an excellent description of the co-added polarized spectra, with a PTE of 12.8% in TE and 34.6% in EE. This already suggests that extensions with, e.g., isocurvature modes can be severely constrained.

Despite this overall agreement, we note that the PTEs are not uniformly good for all cross-frequency spectra (see in particular the 100 × 100 and 100 × 217 in TE). This shows that the baseline instrumental model needs to include further effects to describe all of the data in detail, even if the averages over frequencies appear less affected. The green line in Fig. 28 (mostly visible in the ΔCEE\hbox{$\Delta C_\ell^{EE}$} plot) shows the best-fit leakage correction (shown on its own in Fig. 23), which is obtained when fixing the cosmology to the TT-based model. Let us recall, though, that this correction is for illustrative purposes only, and it is set to zero for all actual parameter searches. Indeed, we shall see that these leakage effects are not enough to bring all the data into full concordance with the model.

In more quantitative detail, Fig. 30 shows the binned (Δ = 100) residuals for the co-added CMB spectra in units of the standard deviation of each data point, (data  model)/error. For TT, we find the greatest deviations at ≈ 434 (−1.8σ), 464(2.7σ), 1214 (−2.1σ), and 1450 (−1.8). At = 1754, where we previously reported a deficit due to the imperfect removal of the 4He-JT cooler line (see Planck Collaboration XIII 2016, Sect. 3), there is now a less significant fluctuation, at the level of −1.4σ. The residuals in polarization show similar levels of discrepancy.

In order to assess whether these deviations are specific to one particular frequency channel or appear as a common signal in all the spectra, Fig. 31 shows foreground-cleaned TT power spectra differences across all frequencies, in units of standard deviations (details on how this is derived can be found in Appendix C.3.2). The agreement between TT spectra is clearly quite good. Figure 32 then shows the residuals per frequency for the TT power spectra with respect to the ΛCDM Plik TT+tauprior best-fit model (see also the zoomed-in residual plots in Fig. C.5). The ≈ 434, 464, and 1214 deviations from the model appear to be common to all frequency channels, with differences between the frequencies smaller than 2σ. However, the deviation at ≈ 1450 is higher at 217 × 217 than in the other channels. In particular, the inter-frequency differences (Fig. 31) between the 217 × 217 power spectrum and the 100 × 100, 143 × 143, and 143 × 217 ones show deviations at ≈ 1450 at the roughly 1.7, 2.6, and 3.4σ levels, respectively.

This inter-frequency difference is due to a deficit in the residuals of the 217 × 217 channel of about −3.4σ in the bin centred at 1454 in Fig. 32. To better quantify this deviation, we also fit for a feature of the type cos2((π/ 2)(p) / (Δ)), with maximum amplitude centred at p = 1460, width Δ = 25 (we impose the feature to be zero at | p | > Δ) and with an independent amplitude in each frequency channel. At 217 × 217, we find an amplitude of (− 37.44 ± 9.5)μK2, while in the other channels we find (− 15.0 ± 7.8)μK2 at 143 × 143 and (− 19.7 ± 7.9)μK2 at 143 × 217. This outlier seems to be at least in part due to chance correlation between the CMB and dust. Indeed, the amplitude of the feature in the different spectra is in rough agreement with the dust emission law. Moreover, the feature can also be found when varying the retained sky fraction in the galactic mask, again with an amplitude scaling compatible with a dust origin. We discuss below the impact on cosmological parameters, see the case “CUT = 1404−1504” in Fig. 35.

Finally we note that there is a deficit in the = 500−800 region (in particular between = 700 and 800) in the residuals of all the frequency spectra, roughly in correspondence with the position of the second and third peaks. Section 4.1 is dedicated to the study of these deviations and their impact on cosmological parameters. In spite of these marginally significant deviations from the model, the χ2 values shown in Table 16 indicate that the ΛCDM model is an acceptable fit to each of the unbinned individual frequency power spectra, with PTEs always \hbox{$\mathcal{P}\gtrsim 10\%$} in TT. We therefore proceed to examine the parameters of the best-fit model.

thumbnail Fig. 32

Residuals in the half-mission TT power spectra after subtracting the Plik TT+tauprior ΛCDM best-fit model (blue points, except for those which differ by at least 2 or 3σ, which are coloured in orange or red, respectively). The light blue line shows the difference between the best-fit model obtained assuming a ΛCDM+AL model and the ΛCDM best-fit baseline; the green line shows the difference of best-fit models using the max = 999 likelihood (fixing the foregrounds to the baseline solution) minus the baseline best-fit (both in the ΛCDM framework); while the pink line is the same as the green one but for max = 1404 instead of max = 999; see text in Sect. 4.1. For the TE and EE spectra, see Fig. 40.

thumbnail Fig. 33

ΛCDM parameters posterior distribution for Plik TT+tauprior. The lower left triangle of the matrix displays how the constraints are modified when the information from one of the frequency channels is dropped. The upper right triangle displays how the constraints are modified when the information from multipoles greater or less than 1000 is dropped. All the results shown in this figure were obtained using the CAMB code.

thumbnail Fig. 34

TE (left) and EE (right) residuals conditioned on the TT spectrum (black line) with 1 and 2σ error bands. The blue points are the actual TE and EE residuals. We do not include any beam-leakage correction here.

The cosmological parameters of interest are summarized in Table 17. Let us note that the cosmological parameters inferred here are obtained using the same codes, priors, and assumptions as in Planck Collaboration XIII (2016), except for the fact that we use the much faster PICO (Fendt & Wandelt 2007a) code instead of CAMB when estimating cosmological parameters12 from TT,TE or TT, TE, EE using high-Planck data. Appendix C.5 establishes that the results obtained with the two codes only differ by small fractions of a standard deviation (less than 15% for most parameters, with a few more extreme deviations). However, we still use the CAMB code for results from EE alone, since in this case the parameter space explored is so wide that it includes regions outside the PICO interpolation region (see Appendix C.5 for further details).

Figure 33 shows the posterior distributions of each pair of parameters of the base ΛCDM model from Plik TT+tauprior. The upper-right triangle compares the 1σ and 2σ contours for the full likelihood with those derived from only the < 1000 or the ≥ 1000 data. Section 4.1.6 addresses the question of whether the results from these different cases are consistent with what can be expected statistically. The lower-left triangle further shows that the results are not driven by the data from a specific channel, i.e., dropping any of the 100, 143, or 217 GHz map data from the analysis does not lead to much change. The next section provides a quantitative analysis of this and other jack-knife tests.

Table 16

Goodness-of-fit tests for the Plik temperature and polarization spectra at high .

Table 17

Cosmological parameters used in this analysis.

We now turn to polarization results. Inter-frequency comparisons and residuals for TE and EE spectra are analysed in detail in Sect. 4.4. Suffice it to say here that the results are less satisfactory than in TT, both in the consistency between frequency spectra and in the detailed χ2 results. This shows that the instrumental data model for polarization is less complete than for temperature, with residual effects at the μK2 level. The model thus needs to be further developed to take full advantage of the HFI data in polarization, given the level of noise achieved. We thus consider the high- polarized likelihood as a “beta” version. Despite these limitations, we include it in the product delivery, to allow external reproduction of the results, even though the tests that we show indicate that it should not be used when searching for weak deviations (at the μK2 level) from the baseline model.

Nevertheless, we generally find agreement between the TT, TE, and EE spectra. Figure 34 shows the TE, and EE residual spectra conditioned on TT, which are close to zero. This is particularly the case for TE below = 1000, which gives some confidence in the polarization model. Most of the data points for TE and EE lie in the ± 2σ range. As for all χ2-based evaluations, the interpretation of this result depends crucially on the quality of the error estimates, i.e., on the quality of our noise model (see Sect. 3.4.4). We further note that the agreement is consistent with the finding that unmodelled instrumental effects in polarization are at the μK2 level.

4. Assessment of the high-multipole likelihood

This section describes tests that we performed to assess the accuracy and robustness of the reference results of the high- likelihood that were presented above. First we establish the robustness of the TT results using Plik alone in Sect. 4.1 and with other likelihoods in Sect. 4.2. We verify in Sect. 4.3 that the amplitudes of the compact-source contributions derived at various frequencies are consistent with our current knowledge of source counts. We then summarize in Sect. 4.4 the results of the detailed tests of the robustness of the polarization results, which are expanded upon in Appendix C.3.5. The paper Planck Collaboration XVI (2016) examines the dependence of the power spectrum on angular direction.

4.1. TT robustness tests

Figure 35 shows the marginal mean and the 68% CL error bars for cosmological parameters calculated assuming different data choices, likelihoods, parameter combinations, and data combinations. The 31 cases shown assume a base-ΛCDM framework, except when otherwise specified. The reference case uses the Plik TT+tauprior data combination. Figure 36 adds the specific results for the lensing parameter AL (left) in a ΛCDM+ AL framework and for the effective number of relativistic species Neff (right) in a ΛCDM+ Neff extended framework.

In both figures, the grey bands show the standard deviation of the parameter shifts relative to the baseline likelihood expected when using a sub-sample of the data (e.g., excising -ranges or frequencies). Because the data sets used to make inferences about a model are changed, one would naturally expect the inferences themselves to change, simply because of the effects of noise and cosmic variance. The inferences could also be influenced by inadequacies in the model, deficiencies in the likelihood estimate, and systematic effects in the data. Indeed, one may compare posterior distributions from different data subsets with each other and with those from the full data set, in order to assess the overall plausibility of the analysis.

To this end it is useful to have some idea about the typical variation in posteriors that one would expect to see even in the ideal case of an appropriate model being used to fit data sets with correct likelihoods and no systematic errors. It can be shown (Gratton & Challinor, in prep.) that if Y is a subset of a data set X, and PX and PY are vectors of the maximum-likelihood parameter values for the two data sets, then the sampling distribution of the differences of the parameter values is given by (PYPX)(PYPX)T=cov(PY)cov(PX),\begin{equation} \overline{\left(\vec{P}_{Y}-\vec{P}_{X}\right)\left(\vec{P}_Y-\vec{P}_{X}\right)^\tens{T}}= \mathrm{cov}(\vec{P}_Y)- \mathrm{cov}(\vec{P}_X), \label{eq:covresult} \end{equation}(53)i.e., the covariance of the differences is simply the difference of their covariances. Here the covariances are approximated by the inverses of the appropriate Fisher information matrices evaluated for the true model. One might thus expect the scatter in the modes of the posteriors to follow similarly, and to be able, if the parameters are well-constrained by the data, to use covariances of the appropriate posteriors on the right-hand side.

thumbnail Fig. 35

Marginal mean and 68% CL error bars on cosmological parameters estimated with different data choices for the Plik likelihood, in comparison with results from alternate approaches or model. We assume a ΛCDM model and use variations of the Plik TT likelihood in most of the cases, in combination with a prior τ = 0.07 ± 0.02 (using neither low- temperature nor polarization data). The “Plik TT+tauprior” case (black dot and thin horizontal black line) indicates the baseline (HM, min = 30, max = 2508), while the other cases are described in Sect. 4.1 (and 4.2, 5.6, E.4). The grey bands show the standard deviation of the expected parameter shift, for those cases where the data used is a sub-sample of the baseline likelihood (see Eq. (53)). All the results were run with PICO except for few ones that were run with CAMB , as indicated in the labels.

4.1.1. Detset likelihood

We have verified (case “DS”) that the results obtained using the half-mission cross-spectra likelihood are in agreement with those obtained using the detset (DS) cross-spectra likelihood. As explained in Sect. 3.4.4, the main difficulty in using the DS likelihood is that the results might depend on the accuracy of the correlated noise correction. Reassuringly, we find that the results from the HM and DS likelihoods agree within 0.2σ. This is an important cross-check, since we expect the two likelihoods to be sensitive to different kinds of temporal systematics. Direct differences of half-mission versus detset-based TT cross-frequency spectra are compared in Fig. 31 (Fig. 41 shows similar plots for the TE and EE spectra.).

When using the detsets, we fit the calibration coefficients of the various detector sets with respect to a reference. The resulting best-fit values are very close to one13, with the greatest calibration refinement being less than 0.2%, in line with the accuracy expected from the description of the data processing in Planck Collaboration VIII (2016). This verifies that the maps produced by the HFI DPC and used for the half-mission-based likelihood come from the aggregation of well-calibrated and consistent data.

4.1.2. Impact of Galactic mask and dust modelling

We have tested the robustness of our results with respect to our model of the Galactic dust contribution in various ways.

Galactic masks.

We have examined the impact of retaining a smaller fraction of the sky, less contaminated by Galactic emission. The baseline TT likelihood uses the G70, G60, and G50 masks (see Appendix A) at 100, 143, and 217 GHz, respectively. We have tested the effects of using G50, G41, and G41 (corresponding to fskynoap=0.60\hbox{$\fsky^\textrm{noap}=0.60$}, 0.50, and 0.50 before apodization, case “M605050” in Fig. 35), and of the priors on the Galactic dust amplitudes relative to these masks described in Table 11. We find stable results as we vary these sky cuts, with the greatest shift in θMC of 0.5σ, compatible with the expected shift of 0.57σ calculated using Eq. (53). Going to higher sky fraction is more difficult. Indeed, the improvement in the parameter determination from increasing the sky fraction at 143 GHz and 217 GHz would be modest, as we would only gain information in the small-scale regime, which is not probed by 100 GHz. Increasing the sky fraction at 100  GHz is also more difficult because our estimates have shown that adding as little as 5% of the sky closer to the Galactic plane requires a change in the dust template and more than doubles the dust contamination at 100  GHz.

Amplitude priors.

We have tested the impact of not using any prior (i.e., using arbitrarily wide, uniform priors) on the Galactic dust amplitudes (case “No gal. priors” in Fig. 35). Again, cosmological results are stable, with the greatest shifts in ln(1010As) of 0.23σ and in ns of 0.20σ. The values of the dust amplitude parameters, however, do change, and their best-fit values increase by about 15 μK2 for all pairs of frequencies, while at the same time the error bars of the dust amplitude parameters increase very significantly. All of the amplitude levels obtained from the 545 GHz cross-correlation are within 1σ of this result. The dust levels from this experiment are clearly unphysically high, requiring \hbox{$22\,\muKsq\ (\mathcal{D}_\ell, \ell=200)$} for the 100 × 100 pair. This level of dust contamination is clearly not allowed by the 545 × 100 cross-correlation, demonstrating that the prior deduced from it is informative. Nevertheless, the fact that cosmological parameters are barely modified in this test indicates that the values of the dust amplitudes are only weakly correlated with those of the cosmological parameters, consistent with the results of Figs. 44 and 45 below, which show the parameter correlations quantitatively.

Galactic dust template slope.

We have allowed for a variation of the Galactic dust index n, defined in Eq. (24), from its default value n = −2.63, imposing a Gaussian prior of −2.63 ± 0.05 (“GALINDEX” case in Fig. 35). We find no shift in cosmological parameters (smaller than ~ 0.1σ) and recover a value for the index of n = −2.572 ± 0.038, consistent with our default choice.

Impact of 500 at 217 GHz.

We have analysed the impact of excising the first 500 multipoles (“LMIN=505 at 217 GHz” in Fig. 35) in the 143 × 217 and 217 × 217 spectra, where the Galactic dust contamination is the strongest. We find very good stability in the cosmological parameters, with the greatest change being a 0.16σ increase in ns. This is compatible with the expectations estimated from Eq. (53) of 0.14σ. The inclusion of the first 500 multipoles at 217 GHz in the baseline Plik likelihood is one of the sources of the roughly 0.45σ difference in ns observed when using the CamSpec code, since the latter excises that range of multipoles; for further discussion see Planck Collaboration XIII (2016, Table 1 and Sect. 3.1), as well as Sect. 4.2.

thumbnail Fig. 36

Marginal mean and 68% CL error bars on the parameters AL (left) and Neff (right) in ΛCDM extensions, estimated with different data choices for the Plik TT likelihood in comparison with results from alternate approaches or model, combined with a Gaussian prior on τ = 0.07 ± 0.02 (i.e., neither low- temperature nor polarization data). The “Plik TT+tauprior” case indicates the baseline (HM, min = 30, max = 2508), while the other cases are described in subsections of Sect. 4.1. The thin horizontal black line shows the baseline result and the thick dashed grey line displays the ΛCDM value (AL = 1 and Neff = 3.04). The grey bands show the standard deviation of the expected parameter shift, for those cases where the data used is a sub-sample of the baseline likelihood (see Eq. (53)).

4.1.3. Impact of beam uncertainties

The case labelled “BEIG” in Fig. 35 corresponds to the exploration of beam eigenvalues with priors 10 times higher than indicated by the analysis of our MC simulation of beam uncertainties (which indicated by dotted lines in Fig. 21). This demonstrates that these beam uncertainties are so small in this data release that they do not contribute to the parameter posterior widths. They are therefore not enabled by default.

4.1.4. Inter-frequency consistency and redundancy

We have tested the effect of estimating parameters while excluding one frequency channel at a time. In Figs. 33 and 35, the “no100” case shows the effect of excluding the 100 × 100 frequency spectrum, the “no143” of excluding the 143 × 143 and 143 × 217 spectra, and the “no217” of excluding the 143 × 217 and 217 × 217 spectra.

We obtain the greatest deviations in the “no217” case for ln(1010As) and τ, which shift to lower values by 0.53σ and 0.47σ, about twice the expected shift calculated using Eq. (53), 0.25σ and 0.23σ respectively (in units of standard deviations of the “no217” case). The value of Ωch2 decreases by only −0.1σ. Figure 37 further shows the 217 × 217 spectrum conditioned on the 100 × 100 and 143 × 143 ones. This conditional deviates significantly in two places, at = 200 and = 1450. The = 1450 case was already discussed in Sect. 3.8 and is further analysed in Sect. 4.1.6. Around = 200, we see some excess scatter (both positive and negative) in the data around a jump between two consecutive bins of the conditional. This corresponds to the two bins around the first peak (one right before and the other almost at the location of the first peak), as can be seen in Fig. 28. All of the frequencies exhibit a similar behaviour (see Fig. 32); however, it is most pronounced in the 217 GHz case. This multipole region is also near the location of the bump in the effective dust model. The magnitude of this excess power in the model is not big enough or sharp enough to explain this excess scatter (see Fig. 17). Finally, note that the best-fit CMB solution at large scales is dominated by the 100 × 100 data, which are measured on a greater sky fraction (see Fig. 14).

This test shows that the parameters of the ΛCDM model do not rely on any specific frequency map, except for a weak pull of the higher resolution 217 GHz data towards higher values of both As and τ (but keeping Asexp(−2τ) almost constant).

thumbnail Fig. 37

217 × 217 spectrum conditioned on the joint result from the 100 × 100 and 143 × 143 spectra. The most extreme outliers are at = 200 and = 1450.

4.1.5. Changes of parameters with min

We have checked the stability of the results when changing min from the baseline value of min = 30 to min = 50 and 100 (and min = 1000, which is discussed in Sect. 4.1.6). These correspond to the cases labelled “LMIN 50” and “LMIN 100” in Fig. 35 (to be compared to the reference case “Plik TT+tauprior”). This check is important, since the Gaussian approximation assumed in the likelihood is bound to fail at very low (for further discussion, see Sect. 3.6).

The results are in good agreement, with shifts in parameters smaller than 0.2σ, well within expectations calculated from Eq. (53). This is also confirmed in Fig. 42, where the TT hybridization scale of the full likelihood is varied (i.e., the multipole where the low- and high- likelihoods are joined).

4.1.6. Changes of parameters with max

We have tested the stability of our results against changes in the maximum multipole max considered in the analysis. We test the restriction to max in the range max = 999−2310, with the baseline likelihood having max = 2508. For each frequency power spectrum we choose maxfreq=min(max,maxfreq,base)\hbox{$\lmax ^{\mathrm{freq}}=\mathrm{min}(\lmax ,\lmax ^{\mathrm{ freq,\,base} })$}, where maxfreq,base\hbox{$\lmax ^{\mathrm{freq,\, base}} $} is the baseline max at each frequency as reported in Table 16. The results shown in Fig. 35 use the same settings as the baseline likelihood (in particular, we leave the same nuisance parameters free to vary) and always use a prior on τ.

The results in Fig. 35 suggest there is a shift in the mean values of the parameters when using low max; e.g., for max = 999, ln(1010As), τ, and Ωch2 are lower by 1.0, 0.8, and 0.8σ with respect to the baseline parameters. These parameters then converge to the baseline values for max ≳ 1500. Following the arguments given earlier (Eq. (53)), when using these nested sub-samples of the baseline data we expect shifts of the order of 0.5, 0.4, and 0.8σ respectively, in units of the standard deviation of the max = 999 results. We further note that the value of θ for max ≲ 1197 is lower compared to the baseline value. In particular, at max = 1197, its value is 0.8σ low, while the expected shift is of the order of 0.7σ, in units of the standard deviation of the max = 1197 results. The value of θ then rapidly converges to the baseline for max ≳ 1300. Figure C.8 in Appendix C.3.3 also shows that these shifts are related to a change in the amplitude of the foreground parameters. In particular, the overall level of foregrounds at each frequency decreases with increasing max, partially compensating for the increase in ln(1010As) and Ωch2. Although all these shifts are compatible with expectations within a factor of 2, we performed some further investigations in order to understand the origin of these changes. In the following, we provide a tentative explanation.

Table 18

Difference of χ2 values between pairs of best-fit models in different ranges for the co-added TT power spectrum.

Table 18 shows the difference in χ2 between the best-fit model obtained using max = 999 (or max = 1404) and the baseline Plik TT+tauprior best-fit solution in different multipole intervals. For this test, we ran the max cases fixing the nuisance parameters to the baseline best-fit solution. This is required in order to be able to “predict” the power spectra at multipoles higher than max, since otherwise the foreground parameters, which are only weakly constrained by the low- likelihood, can converge to unreasonable values. We note that fixing the foregrounds has an impact on cosmological parameters, which can differ from the ones shown in Fig. 35 (see Appendix C.3.4 for a direct comparison). Nevertheless, since the overall behaviour with max is similar, we use this simplified scenario to study the origin of the shifts.

The χ2 differences in Table 18 indicate that the cosmology obtained using max = 999 is a better fit in the region between = 630 and 829. In particular, the low value of θ preferred by the max = 999 data set shifts the position of the third peak to smaller scales. This enables a better fit to the low points at ≈ 700−850 (before the third peak), followed by the high points at ≈ 850−950 (after the third peak). This is also clear from the residuals and the green solid line in Fig. 32, which shows the difference in best-fit models between the max = 999 case and the reference case. However, the values in Table 18 also show that the max = 999 cosmology is disfavoured by the multipole region between ≈ 1330−1430, before the fifth peak. The max = 999 model predicts too little power in this multipole range, which can be better fit if the position of the fifth peak moves to lower multipoles. As a consequence, θ shifts to higher values when including max ≳ 1400.

Concerning the shifts in Ωch2, As and τ, Fig. 35 shows that these parameters converge to the full baseline solution between max = 1404 and max = 1505. The Δχ2 values in Table 18 between the best-fit max = 1404 case and the baseline suggest that the max = 1404 cosmology is disfavoured by the multipole region = 1430−1530 (fifth peak), and – at somewhat lower significance – by the regions close to the fourth peak ( ≈ 1130−1230) and the sixth peak ( ≈ 1730−1829). The pink line in Fig. 32 shows the differences between the max = 1404 best-fit model and the baseline, and it suggests that the max = 1404 cosmology predicts an amplitude of the high- peaks that is too large.

This effect can be compensated by more lensing, which can be obtained with greater values of Ωch2 and ln(1010As), as well as a greater value of τ to compensate for the increase in As in the normalization of the spectra, as observed when considering max ≳ 1500. This also explains why the baseline (max = 2508) best-fit solution prefers a value of the optical depth which is 0.8σ higher than the mean value of the Gaussian prior (τ = 0.07 ± 0.02), τ = 0.085 ± 0.018. In order to verify this interpretation, we performed the following test (using the CAMB code instead of PICO ). We fixed the theoretical lensing power spectrum to the best-fit parameters preferred by the max = 1404 cosmology, and estimated cosmological parameters using the baseline likelihood. This is the “CAMB, FIX LENS” case in Fig. 35, which shows that cosmological parameters shift back to the values preferred at max = 1404 (“CAMB, max = 1404”) if they cannot alter the amount of lensing in the model.

Since the ≈ 1400−1500 region is also affected by the deficit at = 1450 (described in Sect. 3.8), we tested whether excising this multipole region from the baseline likelihood (with max = 2508) has an impact on the determination of cosmological parameters. The results in Fig. 35 (case “CUT = 1404−1504”) show that the parameter shifts are at the level of 0.47, −0.29, 0.38, and 0.45σ on Ωbh2, Ωch2, θ, and ns, respectively (0.39, 0.09, 0.24, and 0.29σ expected from Eq. (53)), confirming that this multipole region has some impact on the parameters, although it cannot completely account for the shift between the max ≈ 1400 case and the baseline.

We also estimated cosmological parameters including only multipoles > 1000 (“LMIN 1000” case), and compared them to the “LMAX 999” case14 (see also Appendix C.5). The two-dimensional posterior distributions in Fig. 33 show the complementarity of the information from ≤ 999 and ≥ 1000, with degeneracy directions between pairs of parameters changing in these two multipole regimes. The min = 1000 likelihood sets constraints on the amplitude of the spectra Ase− 2τ and on ns that are almost a factor of 2 weaker than the ones obtained with the baseline likelihood, and somewhat higher than the ones obtained with max = 999. The value of τ is thus more effectively determined by its prior and shifts downward by 0.59σ with respect to the baseline. The value of Ωch2 shifts upward by 1.7σ (cf. 0.8σ expected from Eq. (53)). Whether this change is just due to a statistical fluctuation is still a matter of investigation.

However, since parameter shifts are correlated, we evaluated whether the ensemble of the shifts in all cosmological parameters between the max = 999 and min = 1000 cases are compatible with statistical expectations. In order to do so, we computed the χΔ2\hbox{$\chi_\Delta^2$} statistic of the shift as χΔ2=ijΔiΣij-1Δj,\begin{equation} \chi_\Delta^2 = \sum_{ij} \Delta_i \Sigma^{-1}_{ij} \Delta_j , \end{equation}(54)where Δi is the difference in best-fit value of the ith parameter between the max = 999 and min = 1000 cases and Σ is the covariance matrix of the expected shifts, calculated as the sum of the parameter covariance matrices obtained in each of the two cases, ignoring correlations between the two datasets. We include in this calculation the ΛCDM parameters (Ωbh2, Ωch2, θ, ns, Asexp(−2τ)), excluding τ, since the constraints on this parameter are dominated by the same prior in both cases, and using Asexp(−2τ) instead of ln(1010As), since the latter is very correlated with τ and the TT power spectrum is mostly sensitive to the combination Asexp(−2τ). Finally, we estimate the χΔ2\hbox{$\chi_\Delta^2$} both in the case where we leave the foregrounds free to vary or in the case where we fix them to the best fit of the baseline PlikTT + tauprior solution. Assuming that χΔ2\hbox{$\chi_\Delta^2$} has a χ2 distribution for 5 degrees of freedom, we find that the shifts observed in the data are consistent with simulations at the 1.2σ (1.1σ with fixed foregrounds) level for the case where we do not include the low-TT likelihood at < 30 to the max = 999 case, and at the 1.5σ (1.4σ with fixed foregrounds) level for the case where we include the low-TT likelihood. We also find that the use of Asexp(−2τ) instead of ln(1010As) changes these significances only in the case where we include the low-TT likelihood to the max = 999 case and leave the foregrounds free to vary, in which case we find consistency at the level of 1.8σ, in agreement with the findings of Addison et al. (2016; although in this case the use of ln(1010As) and the exclusion of τ makes this test less indicative of the true significance of the shifts). In all cases, we do not find evidence for a discrepancy between the two datasets. A more precise and extended evaluation and discussions of these shifts, based on numerical simulations, will be presented in a future publication.

4.1.7. Impact of varying AL

Figure 36 (left) displays the impact of various choices on the value of the lensing parameter AL in the ΛCDM+AL framework. The baseline likelihood prefers a value of AL that is about 2σ greater than the physical value, AL = 1. It is clear that this preference only arises when data with max ≳ 1400 are included, and it is caused by the same effects as we proposed in Sect. 4.1.6 to explain the shifts in parameters at max ≳ 1400 in the ΛCDM case. More lensing helps to fit the data in the ≈ 1300−1500 region, as indicated by the χ2 differences between the ΛCDM+AL best-fit and the ΛCDM one in Table 18. This drives the value of AL to 1.159 ± 0.090 with Plik TT+tauprior, 1.8σ higher than expected. The case “ΛCDM+AL” of Fig. 35 also shows that opening up this unphysical degree of freedom shifts the other cosmological parameters at the 1σ level; e.g., Ωch2 and As shift closer to the values preferred in the ΛCDM case when using max ≲ 1400. While in the ΛCDM case high values of these parameters allow increasing lensing, in the ΛCDM+AL case this is already ensured by a high value of AL, so Ωch2 and As can adopt values that better fit the ≲ 1400 range. When using Plik TT in combination with the lowTEB likelihood, the deviation increases to 2.4σ, AL = 1.204 ± 0.086,15 due to the fact that more lensing allows smaller values of Ωch2 and As and a greater value of ns, better fitting the deficit at ≈ 20 in the temperature power spectrum (see Planck Collaboration XIII 2016, Sect. 5.1.2 and Fig. 13).

4.1.8. Impact of varying Neff

We have investigated the effect of opening up the Neff degree of freedom in order to assess the robustness of the constraints on the ΛCDM extensions, which rely heavily on the high- tail of the data. Figure 36 (right) shows that Neff departs from the standard 3.04 value by about 1σ when using Plik TT+tauprior, Neff = 2.7 ± 0.33. The χ2 improvement for this model over ΛCDM is only Δχ2 = 1.5. We note that when the lowTEB likelihood (or alternatively, the low-TT likelihood plus the prior on τ) is used in combination with Plik TT, the value of Neff shifts higher by about 1σ, Neff = 3.09 ± 0.29. This shift is about a factor 2 more than the one expected from Eq. (53), 0.5σ, between the Plik TT+tauprior and Plik TT+tauprior+low-TT cases. This shift is due to the fact that the deficit at ≈ 20 is better fit by higher ns and, as a consequence, an increase in Neff helps decreasing the enhanced power at high .

Figure 36 also shows that, not surprisingly, the most extreme variations as compared to the reference case (less than 1σ) arise when the high-resolution data are dropped (by reducing max or by removing the 217 GHz channel), owing to the strong dependence of the Neff constraints on the damping tail.

Having opened up this degree of freedom, the standard parameters are now about 1σ away (see case “ΛCDM+Neff” of Fig. 35), and such a model would prefer quite a low value of H0, which would then be at odds with priors derived from direct measurements (see Planck Collaboration XIII 2016, for an in-depth analysis).

4.2. Intercomparison of likelihoods

In addition to the baseline high-Plik likelihood, we have developed four other high- codes, CamSpec , Hillipop , Mspec , and Xfaster . CamSpec and Xfaster have been described in separate papers (Planck Collaboration XV 2014; Rocha et al. 2011), and brief descriptions of Mspec and Hillipop are given in Appendix D. These codes have been used to perform data consistency tests, to examine various analysis choices, and to cross-check each other by comparing their results and ensuring that they are the same. In general, we find good agreement between the codes, with only minor differences in cosmological parameters.

The CamSpec , Hillipop , and Mspec codes are, like Plik , based on pseudo-C estimators and an analytic calculation of the covariance (Efstathiou 2004, 2006), with some differences in the approximations used to calculate this covariance. The Xfaster code (Rocha et al. 2011) is an 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), with noise bias estimated using difference maps, as described in Planck Collaboration IX (2016). For temperature, all of the codes use the same Galactic masks, but they differ in point-source masking: Hillipop uses a mask based on a combination of S/N> 7 and cuts based on flux, while the others use the baseline S/N> 5 mask described in Appendix A. The codes also differ in foreground modelling, in the choice of data combinations, and in the -range. For the comparison presented here, all make use of half-mission maps.

thumbnail Fig. 38

Comparison of power spectra residuals from different high- likelihood codes. The figure shows “data/calib FG Plik CMB”, where “data” stands for the empirical cross-frequency spectra, “FG” and “calib” are the best-fit foreground model and recalibration parameter for each individual code at that frequency, and the best-fit model Plik CMB is subtracted for visual presentation. These plots thus show the difference in the amount of power each code attributes to the CMB. The power spectra are binned in bins of width Δ = 100. The y-axis scale changes at = 500 for TT and = 1000 for EE (vertical dashes).

thumbnail Fig. 39

Comparison of error bars from the different high- likelihood codes. The quantities plotted are the ratios of each code’s error bars to those from Plik , and are for bins of width Δ = 100. Results are shown only in the range common to Plik and the code being compared.

Figure 38 shows a comparison of the power spectra residuals and error bars from each code, while Fig. E.5 in Appendix E.4 compares the combined spectra with the best-fit model. In temperature, the main feature visible in these plots is an overall nearly constant shift, up to 10 μK2 in some cases. This represents a real difference in the best-fit power each code attributes to foregrounds. For context, it is useful to note the statistical uncertainty on the foregrounds; for example, the 1σ error on the total foreground power at 217 GHz at = 1500 is 2.5 μK2 (calculated here with Mspec , but similar for the other codes). Shifts of this level do not lead to very large differences in cosmological parameters except in a few cases that we discuss.

For easier visual comparison of error bars, we show in Fig. 39 the ratios of each code’s error bars to those from Plik . These have been binned in bins of width Δ = 100, and are thus sensitive to the correlation structure of each code’s covariance matrix, up to 100 multipoles into the off-diagonal. For all the codes and for both temperature and polarization, the correlation between multipoles separated by more than Δ = 100 is less than 3%, so Fig. 39 contains the majority of the relevant information about each code’s covariance.

A few differences are visible, mostly at high frequency, when the 217 GHz data are used. First, the Hillipop error bars in TT for 143 × 217 become increasingly tighter than the other codes at > 1700. This is because Hillipop , unlike the other codes, gives non-zero weight to 143 × 217 spectra when both the 143 and the 217 GHz maps come from the same half-mission. This leads to a slight increase in power at high compared to Plik , as can be seen in Fig. 38. Conversely, the Hillipop error bars are slightly larger by a few percent at < 1700; however the source of this difference is not understood. Second, the Mspec error bars in temperature are increasingly tighter towards higher frequency, as compared to other codes; for 217 × 217, Mspec uncertainties are smaller by 510% for between 1000 and 2000. This arises from the Mspec map-based Galactic cleaning procedure, which removes excess variance due to CMB–foreground correlations by subtracting a scaled 545 GHz map. However, for polarization, where one must necessarily clean with the noisier 353 GHz maps, the Mspec error bars for TE and EE become larger. CamSpec , which also performs a map cleaning for low- polarization, switches to a power-spectrum cleaning at higher to mitigate this effect.

The differences in ΛCDM parameters from TT are shown in Table 19. Generally, parameters agree to within a fraction ofσ, but with some differences we discuss. One thing to keep in mind in interpreting this comparison is that these differences are not necessarily indicative of systematic errors. Some of the differences are expected due to statistical fluctuations because different codes weight the data differently.

One of the biggest differences with respect to the baseline code is in ns, which is higher by about 0.45σ for CamSpec , with a related downward shift of Ase− 2τ. To put these shifts into perspective, we refer to the whisker plots of Figs. 35 and 36 which compare CamSpec TT results with Plik in the ΛCDM case (base and extended). A difference in ns of about 0.16σ between Plik and CamSpec can be attributed to the inclusion in Plik of the first 500 multipoles for 143 × 217 and 217 × 217; these multipoles are excluded in CamSpec (see also Sect. 4.1.2). Indeed, cutting out those multipoles in Plik brings ns closer by 0.16σ to the CamSpec value and slightly degrades the constraint on ns compared to the full Plik result. Using Eq. (53), we see that the shift and degradation in constraining power are consistent with expectations. A similar 0.16σ shift can be attributed to different dust templates. CamSpec uses a steeper power law index (−2.7). Using the CamSpec template in Plik brings ns closer to the CamSpec value. Allowing the power law index of the galactic template to vary when exploring cosmological parameters yields a slightly shallower slope (see Sect. 4.1.2). The slope of the dust template is mainly determined at relatively high , i.e., in the regime where it is hardest to determine the template accurately since the dust contribution is only a small fraction of the CIB and point-source contributions (see the ≳ 1000 parts of Figs. 19 and 20). The remaining difference of 0.13σ arises from differences in data preparation (maps, calibration, binning) and covariance estimates. We therefore believe that a 0.2σ is a conservative upper bound of the systematic error in ns associated with the uncertainties in the modelling of foregrounds, which is the biggest systematic uncertainty in TT.

A shift that is less well understood is the ≈ 1σ shift in Ase− 2τ between Plik and Hillipop . The preference for a lower amplitude from Hillipop is sourced by the lower power attributed to the CMB, seen in Fig. 38. With τ partially fixed by the prior, this implies lower As and hence a smaller lensing potential envelope, explaining the somewhat lower value of AL found by Hillipop . Tests performed with the same code suggest that 1σ is too great a shift to be explained simply by the different foreground models, so some part of it must be due to the different data weighting; as can be seen in Fig. 39, Hillipop gives less weight to 500 ≲ ≲ 1500, and slightly more outside of this region.

This comparison also shows the stability of the results with respect to the Galactic cleaning procedure. Mspec and Plik use different procedures, yet their parameter estimates agree to better than 0.5σ (see Appendix D.1). But we note that the Plik –CamSpec differences are higher in the polarization case, and can reach 1σ, as can be judged from the whisker plot in polarization of Fig. C.10.

Table 19

Comparison between the parameter estimates from different high- codes.

4.3. Consistency of Poisson amplitudes with source counts

The Poisson component of the foreground model is sourced by shot-noise from astrophysical sources. In this section we discuss the consistency between the measured Poisson amplitudes and other probes and models of the source populations from which they arise. The Poisson amplitude priors that we calculate are not used in the main analysis, because they improve uncertainties on the cosmological parameters by at most 10%, and only for a few extensions; instead they serve as a self-consistency check.

This type of check was also performed in Like13, which we update here by:

  • 1.

    developing a new method for calculating these priors that isaccurate enough to give realistic uncertainties on Poissonpredictions (for the first time);

  • 2.

    including a comparison of more theoretical models;

  • 3.

    taking into account the 2015 point-source masks.

In Like13 the Poisson power predictions were calculated via C=0ScutdSS2dNdS,\begin{eqnarray} \label{eq:clpoissonold} C_\ell = \int_0^{S_{\rm cut}} {\rm d}S \, S^2 \frac{{\rm d}N}{{\rm d}S} , \end{eqnarray}(55)where dN/ dS is the differential number count, Scut is an effective flux-density cut above which sources are masked, and the integral was evaluated independently at each frequency. Although it is adequate for rough consistency checks, Eq. (55) ignores the facts that the 2013 point-source mask was built from a union of sources detected at different frequencies, and that the Planck flux-density cut varies across the sky, and it also ignores the effect of Eddington bias. In order to accurately account for all of these effects, we now calculate the Poisson power as Cij=0dS1...dSnSiSjdN(S1,...,Sn)dS1...dSnI(S1,...,Sn),\begin{eqnarray} \label{eq:clpoisson} C_\ell^{ij} = \int_0^\infty {\rm d}S_1 \ldots {\rm d}S_n \; S_i S_j \frac{{\rm d}N(S_1,\ldots,S_n)}{{\rm d}S_1\ldots {\rm d}S_n} \, I(S_1,\ldots,S_n) , \end{eqnarray}(56)where the frequencies are labelled 1...n, the differential source count model, dN/ dS, is now a function of the flux densities at each frequency, and I(S1,...,Sn) is the joint “incompleteness” of our catalogue for the particular cut that was used to build the point-source mask.

Table 20

Priors on the Poisson amplitudes given a number of different point-source masks and models.

The joint incompleteness was determined by injecting simulated point sources into the Planck sky maps, using the procedure described in Planck Collaboration XXVI (2016). The same point-source detection pipelines that were used to produce the Second Planck Catalogue of Compact Sources (PCCS2) were run on the injected maps, producing an ensemble of simulated Planck sky catalogues with realistic detection characteristics. The joint incompleteness is defined as the probability that a source would not be included in the mask as a function of the source flux density, given the specific masking thresholds being considered. The raw incompleteness is a function of sky location, because the Planck noise varies across the sky. The incompleteness that appears in Eq. (56) is integrated over the region of the sky used in the analysis; the injection pipeline estimates this quantity by injecting sources only into these regions.

Equation (56) can be applied to any theoretical model which makes a prediction for the multi-frequency dN/ dS. We have adopted the following models.

  • 1.

    For radio galaxies we have two models. The first is the Tucciet al. (2011) model, updated toinclude new source-count measurements from Mocanuet al. (2013). We also consider aphenomenological model that is a power law in flux density andfrequency, and assumes that the sources’ spectral indices areGaussian-distributed with mean \hbox{$\bar\alpha$} and standard deviation σα; we use different values for \hbox{$\bar\alpha$} and σα above and below 143 GHz. We shall refer to this second model as the “power-law” model, and the differential source counts are given by dN(S1,S2,S3)dS1dS2dS3=A(S1S2S3)γ12πσ12σ23×exp[(α(S1,S2)α̅12)22σ122(α(S2,S3)α̅32)22σ322],\begin{eqnarray} \label{eq:dNdS} &&\frac{{\rm d}N(S_1,S_2,S_3)}{{\rm d}S_1 {\rm d}S_2 {\rm d}S_3} = \frac{A(S_1 S_2 S_3)^{\gamma-1}}{2\pi \sigma_{12} \sigma_{23}} \\ &&\quad \quad \quad\times \exp \left[-\frac{(\alpha(S_1,S_2) - \bar\alpha_{12})^2}{2\sigma_{12}^2} -\frac{(\alpha(S_2,S_3) - \bar\alpha_{32})^2}{2\sigma_{32}^2} \right] ,\nonumber \end{eqnarray}(57)where labels 13 refer to Planck 100, 143, and 217 GHz and α(Si,Sj) = ln(Sj/Si) / ln(νj/νi). Both radio models are excellent fits to the available source-count data, and we take the difference between them as an estimate of model uncertainty. With the power-law model we are additionally able to propagate uncertainties in the source count data to the final Poisson estimate via MCMC.

  • 2.

    For dusty galaxies we use the Béthermin et al. (2012) model, as in Planck Collaboration XXX (2014). The model is in good agreement with the number counts measured with te Spitzer Space Telescope and the Herschel Space Observatory. It also gives a reasonable CIB redshift distribution, which is important for cross-spectra, and is a very good fit to CIB power spectra (see Béthermin et al. 2013). In contrast to the radio-source case, the major contribution to the dusty galaxy Poisson power arises from sources with flux densities well below the cuts; for example, we note that decreasing the flux-density cuts by a factor of 2 decreases the Poisson power by less than 1% at the relevant frequencies. In this case, Eq. (55) is a sufficient and more convenient approximation, and we make use of it when calculating Poisson levels for dusty galaxies.

We give predictions for Poisson levels for three different masks: (1) the 2013 point-source mask, which was defined for sources detected at S/N> 5 at any frequency between 100 and 353 GHz; (2) the 2015 point-source mask, which is frequency-dependent and includes S/N> 5 sources detected only at each individual frequency (used by Plik , CamSpec , and Mspec in this work); and (3) the Hillipop mask, which is also frequency-dependent and involves both a S/N cut and a flux-density cut16.

Table 20 summarizes the main results of this section. Generally, we find good agreement between the priors from source counts and the posteriors from chains, with the priors being much more constraining. The exception to the good agreement is at 100 GHz where the prediction is lower than the measured value by around 4σ for the baseline 2015 mask and 6σ for the Hillipop mask. This is a sign either of a foreground modelling error or (perhaps more likely) of a residual unmodelled systematic in the data. We note that this disagreement was not present in Like13, where the Poisson amplitude at 100 GHz was found to be smaller. We also note that removing the relative calibration prior (Eq. (35)) or increasing the max at 100 GHz by a few hundred reduces the tension in the Mspec results. In any case, it is unlikely to affect parameter estimates at all, since very little cosmological information comes from the multipole range at 100 GHz that constrains the Poisson amplitude.

4.4. TE and EE test results

4.4.1. Residuals per frequency and inter-frequency differences

Figure 40 shows the residuals for each frequency and Fig. 41 shows the differences between frequencies of the TE and EE power spectra (the procedure is explained in Appendix C.3.2). The residuals are calculated with respect to the best-fit cosmology as preferred by Plik TT+tauprior, although we use the best-fit solution of the Plik TT, TE, EE+tauprior run to subtract the polarized Galactic dust contribution.

thumbnail Fig. 40

Residual frequency power spectra after subtraction of the Plik TT+tauprior best-fit model. We clean Galactic dust from the spectra from using the best-fit solution of Plik TT, TE, EE+tauprior. The residuals are relative to the baseline HM power spectra (blue points, except for those that deviate by at least 2 or 3σ, which are shown in orange or red, respectively). The vertical dashed lines delimit the ranges retained in the likelihood. Upper: TE power spectra. Lower: EE power spectra.

thumbnail Fig. 41

Inter-frequency foreground-cleaned power-spectra differences. Each panel shows the difference of two frequency power spectra, that indicated on the left axis minus that on the bottom axis, after subtracting foregrounds using the best-fit PlanckTT+lowP foreground solutions. Differences are shown for both the HM power spectra (dark blue) and the DS power spectra (light blue).

The binned inter-frequency residuals show deviations at the level of a few μK2 from the best-fit model. These deviations do not necessarily correspond to high values of the χ2 calculated on the unbinned data (see Table 16). This is because some of the deviations are relatively small for the unbinned data and correctly follow the expected χ2 distribution. However, if the deviations are biased (e.g., have the same sign) in some range, they can result in larger deviations (and large χ2) after binning. Thus, the χ2 calculated on unbinned data is not always sufficient to identify these type of biases. We therefore also use a second quantity, χ, defined as the weighted linear sum of residuals, to diagnose biased multipole regions or frequency spectra: χ=wT(C)withw=(diagC)1/2,\begin{eqnarray} \chi = \vec{w}^\tens{T}(\vec{\hat{C}}-\vec{C}) \quad \text{with}\quad \vec{w} = (\mathrm{diag}\,\tens{C})^{-1/2} , \end{eqnarray}(58)where Ĉ is the unbinned vector of data in the multipole region or frequency spectrum of interest, C) is the corresponding model, and w is a vector of weights, equal to the inverse standard deviation evaluated from the diagonal of the corresponding covariance matrix C. The χ statistic is distributed as a Gaussian with zero mean and standard deviation equal to σχ=wTCw.\begin{eqnarray} \sigma_\chi=\sqrt{\vec{w}^\tens{T}\tens{C}\vec{w}} . \end{eqnarray}(59)We then define the normalized χnorm as the χ in units of standard deviation, χnorm=χ/σχ.\begin{equation} \chi_{\mathrm{norm}}=\chi/\sigma_\chi . \label{chinorm} \end{equation}(60)The χnorm values that we obtain for different frequency power spectra are given in Table 16.

For EE, the worst-behaved spectra from the χnorm point of view are 143 × 143 (3.7σ deviation) and 100 × 217 (−3.0σ), while from the χ2 point of view, the worst is 100 × 143 (PTE = 3.9%). For TE, the worst from the χnorm point of view are 100 × 217 (5σ), 100 × 100 (3.7σ), and 143 × 143 (−2.2σ), while from the χ2 point of view the worst is 100 × 100 (PTE = 0.43%). The extreme deviations from the expected distributions show that the frequency spectra are not described very accurately by our data model. This is also clear from Fig. 41, which shows that there are differences of up to 5σ between pairs of foreground-cleaned spectra.

However, as the co-added residuals in Fig. 29 show, systematic effects in the different frequency spectra appear to average out, leaving relatively small residuals with respect to the Plik TT+tauprior best-fit cosmology. In other words, these effects appear not to be dominated by common modes between detector sets or across frequencies. This is also borne out by the good agreement between the data and the expected polarization power spectra conditioned on the temperature ones, as shown in the conditional plots of Fig. 34.

4.4.2. TE and EE robustness tests

For TE and EE, we ran tests of robustness similar to those applied earlier to TT. These are presented in Appendix C.3.5, and the main conclusions are the following. We find that the Plik cosmological results are affected by less than 1σ when using detset cross-spectra instead of half-mission ones. This is also the case when we relax the dust amplitude priors, when we marginalize over beam uncertainties, or when we change min or max. The alternative CamSpec likelihood has larger shifts, but still smaller than 1σ in TE and 0.5σ in EE. However, we also see larger shifts (more than 2σ in TE) with Plik when some frequency channels are dropped; and, when they are varied, the beam leakage parameters adopt much higher values than expected from the prior, while still leaving some small discrepancies between individual cross-spectra that have yet to be explained.

These results shows that our data model leaves residual instrumental systematic errors and is not yet sufficient to take advantage of the full potential of the HFI polarization information. Indeed, the current data model and likelihood code do not account satisfactorily for deviations at the μK2 level, even if they can be captured in part by our beam leakage modelling. Nevertheless, the results for the ΛCDM model obtained from the Plik TE+tauprior and Plik EE+tauprior runs are in good agreement with the results from Plik TT+tauprior (see Appendix C.3.6). This agreement between temperature and polarization results within ΛCDM is not a proof of the accuracy of the co-added polarization spectra and their data model, but rather a check of consistency at the μK2 level. This consistency is, of course, a very interesting result in itself. But this comparison of probes cannot yet be pushed further to check for the potential presence of a physical inconsistency within the base model that the data could in principle detect or constrain.

5. The full Planck spectra and likelihoods

This section discusses the results that are obtained by using the full Planck likelihood. Section 5.1 first addresses the question of robustness with respect to the choice of the hybridization scale (the multipole at which we transition from the low- likelihood to the high- likelihood). Sections 5.2 and 5.3 then present the full results for the power spectra and the baseline cosmological parameters. Section 5.4 summarizes the full systematic error budget. Section 5.5 concentrates on the significance of the possibly anomalous structure around ≈ 20 in this new release. We then introduce in Sect. 5.6 a useful compressed Planck high- temperature and polarization CMB-only likelihood, Plik_lite , which, when applicable, enables faster parameter exploration. Finally, in Sect. 5.7, we compare the Planck 2015 results with the previous results from WMAP, ACT, and SPT.

5.1. Insensitivity to hybridization scale

thumbnail Fig. 42

Marginal mean and 68% CL error bars on cosmological parameters estimated with different multipoles for the transition between the low- and the high- likelihood. Here we use only the TT power spectra and a Gaussian prior on the optical depth τ = 0.07 ± 0.02, within the base-ΛCDM model. “Plik TT+tauprior” refers to the case where we use the Plik high- likelihood only.

Before we use the low- and high- likelihoods together, we address the question of the hybridization scale, hyb, at which we switch from one to the other (neglecting correlations between the two regimes, as we did and checked in Like13). To that end, we focus on the TT case and use a likelihood based on the Blackwell-Rao estimator and the Commander algorithm (Chu et al. 2005; Rudjord et al. 2009) as described in Sect. 2.2, since this likelihood can be used to much higher max than the full pixel-based T,E,B one. For this test without polarization data, we assume the same τ = 0.07 ± 0.02 prior as before.

The whisker plot of Fig. 42 shows the marginal mean and the 68% CL error bars for base-ΛCDM cosmological parameters when hyb is varied from the baseline value of 30 (case “LOWL 30”) to hyb =50, 100, 150, 200, and 250, and compared to the Plik TT+tauprior case. The difference between the “LOWL 30” and “Plik TT+tauprior” values shows the effect of the low- dip at ≈ 20, which reaches 0.5σ on ns. The plot shows that the effect of varying hyb from 30 to 150 is a shift in ns by less than 0.1σ. This is the result of the Gaussian approximation pushed to min = 30, already discussed in the simulation section (Sect. 3.6). It would have been much too slow to run the full low-TEB likelihood with max substantially greater than 30, and we decided against the only other option, to leave a gap in polarization between = 30 and the hybridization scale chosen in TT.

5.2. The Planck 2015 CMB spectra

The visual appearance of Planck 2015 CMB co-added spectra in TT, TE, and EE can be seen in Fig. 50. Goodness-of-fit values can be found in Table E.1 of Appendix E. These differ somewhat from those given previously in Table 16 for Plik alone, because the inclusion of low in temperature brings in the ≈ 20 feature (see Sect. 5.5). Still, they remain acceptable, with PTEs all above 10% (16.8% for TT).

With this release, Planck now detects 36 extrema in total, consisting of 19 peaks and 17 troughs. Numerical values for the positions and amplitudes of these extrema may be found in Table E.2 of Appendix E.2, which also provides details of the steps taken to derive them. We provide in Appendix E.3 an alternate display of the correlation between temperature and (E-mode) polarization by showing their Pearson correlation coefficient and their decorrelation angle versus scale (Figs. E.2 and E.3).

5.3. Planck 2015 model parameters

Figure 43 compares constraints on pairs of parameters as well as their individual marginals for the base-ΛCDM model. The grey contours and lines correspond to the results of the 2013 release (Like13), which was based on TT and WMAP polarization at low (denoted by WP), using only the data from the nominal mission. The blue contours and lines are derived from the 2015 baseline likelihood, Plik TT+lowTEB (“PlanckTT+lowP” in the plot), while the red contours and line are obtained from the full Plik TT, EE, TE+lowTEB likelihood (“PlanckTT, TE, EE+lowP” in the plot, see Appendix E.1 for the relevant robustness tests). In most cases the 2015 constraints are in quite good agreement with the earlier constraints, with the exception of the normalization As, which is higher by about 2%, reflecting the 2015 correction of the Planck calibration which was indeed revised upward by about 2% in power. The figure also illustrates the consistency and further tightening of the parameter constraints brought by adding the E-mode polarization at high . The numerical values of the Planck 2015 cosmological parameters for base ΛCDM are given in Table 21.

thumbnail Fig. 43

ΛCDM parameter constraints. The grey contours show the 2013 constraints, which can be compared with the current ones, using either TT only at high (red) or the full likelihood (blue). Apart from further tightening, the main difference is in the amplitude, As, due to the overall calibration shift.

As shown in Fig. 44, the degeneracies between foreground and calibration parameters generally do not affect the determination of the cosmological parameters. In the Plik TT+lowTEB case (top panel), the dust amplitudes appear to be nearly uncorrelated with the basic ΛCDM parameters. Similarly, the 100 and 217 GHz channel calibration is only relevant for the level of foreground emission. Cosmological parameters are, however, mildly correlated with the point-source and kinetic SZ amplitudes. Correlations are strongest (up to 30%) for the baryon density (Ωbh2) and spectral index (ns). We do not show correlations with the Planck calibration parameter (yP), which is uncorrelated with all the other parameters except the amplitude of scalar fluctuations (As). The bottom panel shows the correlation for the Plik TE+lowTEB and Plik EE+lowTEB cases, which do not affect the cosmological parameters, except for 20% correlations in EE between the spectral index (ns) and the dust contamination amplitude in the 100 and 143 GHz maps.

Table 21

Constraints on the basic six-parameter ΛCDM model using Planck angular power spectra.

We also display in Fig. 45 the correlations between the foreground parameters and the cosmological parameters in the Plik TT+lowTEB case when exploring classical extensions to the ΛCDM model. While nrun seems reasonably insensitive to the foreground parameters, some extensions do exhibit a noticeable correlation, up to 40% in the case of YHe and the point-source level at 143 GHz.

Finally, we note that power spectra and parameters derived from CMB maps obtained by the component-separation methods described in Planck Collaboration IX (2016) are generally consistent with those obtained here, at least when restricted to the < 2000 range in TT; this is detailed in Sect. E.4.

thumbnail Fig. 44

Parameter correlations for Plik TT+lowTEB (top), Plik TE+lowTEB (bottom left), and Plik EE+lowTEB (bottom right). The degeneracies between foreground and calibration parameters do not strongly affect the determination of the cosmological parameters. In these figures the lower triangle gives the numerical values of the correlations in percent (with values below 10% printed at the smallest size), while the upper triangle represents the same values using a colour scale.

thumbnail Fig. 45

Parameter correlations for Plik TT+lowTEB, including some ΛCDM extensions. The leftmost column is identical to Fig. 44 and is repeated here to ease comparison. Including extensions to the ΛCDM model changes the correlations between the cosmological parameters, sometimes dramatically, as can be seen in the case of AL. There is no correlation between the cosmological parameters (including the extensions) and the dust amplitude parameters. In most cases, the extensions are correlated with the remaining foreground parameters (and in particular with the point-source amplitudes at 100 and 143 GHz, and with the level of CIB fluctuations) with a strength similar to those of the other cosmological parameters (i.e., less than 30%). YHe exhibits a stronger sensitivity to the point-source levels.

5.4. Overall systematic error budget assessment

The tests presented throughout this paper and its appendices documented our numerous tests of the Planck likelihood code and its outputs. Here, we summarize those results and attempt to isolate the dominant sources of systematic uncertainty. This assessment is of course a difficult task. Indeed, all known systematics are normally corrected for, and when relevant, the uncertainty on the correction is included in the error budget and thus in the error bar we report. In that sense, except for a very few cases where we decided to leave a known uncertainty in the data, this section tries to deal with the more difficult task of evaluating the unknown uncertainty!

This section summarizes the contribution of the known systematic uncertainties along with these potential unknown unknowns, specifically highlighting both internal consistency tests based on comparing subsets of the data, along with those using end-to-end instrumental simulations.

5.4.1. Low- budget

The low- likelihood has been validated using both internal consistency tests and simulation-based, tests. Here we summarize only the main result of the analysis, which has been set forth in Sect. 2 above.

A powerful consistency test of the polarization data, described in Sect. 2.4, is derived by rotating some of the likelihood components by π/ 4. Specifically, the rotation is applied to the data maps and only to the noise covariance matrix (the likelihood being a scalar function, applying the same rotation to the signal matrix as well would be equivalent to not performing the rotation). The net effect is a conversion of E → −B and BE for the signal, but leaving unaffected the Gaussian noise described in the covariance matrix. Under these circumstances, we do not expect to pick up any reionization signal, since it would then be present in BB or TB: the operation should result in a null τ detection. This is precisely what happens (see the blue dashed curve in Fig. 8). It is of course possible – though unlikely – that systematics are only showing up in the E channel, leaving B modes unaffected. Indeed, this possibility is further challenged by the fact that we do not detect anomalies in any of the six polarized power spectra; as detailed in Fig. 7, they are consistent with a ΛCDM signal and noise as described by the final 70 GHz covariance matrix.

These tests are specific to Planck and aimed at validating the internal consistency of the datasets employed to build the likelihood. As a further measure of consistency, we have carried out a null test employing the WMAP data, detailed in Sect. 2.6. In brief, we have taken WMAP’s Ka, Q and V channels and cleaned them from any polarized foreground contributions using a technique analogous to the one used to clean the LFI 70 GHz maps, employing the Planck 353 GHz map to minimize any dust contribution, but relying on WMAP’s K channel to remove any synchrotron contribution. The resulting LFI 70 GHz and WMAP maps separately lead to compatible τ detections; their half-difference noise estimates are compatible with their combined noise and do not exhibit a reionization signal, as shown in Fig. 11.

We learn from these tests that if the EE and TE signal we measure at 70 GHz is due to systematics, then these systematics should affect only the above spectra in such a way to mimic a genuine reionization signal, and one that is fully compatible in the maps with that present in (cleaned) WMAP data. This is extremely unlikely and conclude that Planck 70 GHz is dominated by a genuine contribution from the sky, compatible with a signal from cosmic reionization.

The tests described so far do not let us accurately quantify the magnitude of a possible systematic contribution, nor to exclude artefacts arising from the data pipeline itself and, specifically, from the foreground cleaning procedure. These can be only controlled through detailed end-to-end tests, using the FFP8 simulations (Planck Collaboration XII 2016). As detailed in Sect. 2.5, we have performed end-to-end validation with 1000 simulated frequency maps containing signal, noise, and foreground contributions as well as specific systematics effects, mimicking all the steps in the actual data pipeline. Propagating to cosmological parameters (τ and As, which are most relevant at low in the ΛCDM model), we detect no bias within the simulation error budget. The total impact of any unknown systematics on the final τ estimate is at most 0.1σ. This effectively rules out any detectable systematic contribution from the data pipeline or or from the instrumental effects considered in the FFP8 simulations. A complementary analysis has been performed in Planck Collaboration III (2016), including further systematic contributions not incorporated into FFP8. This study, which should be taken as a worst-case scenario, limits the possible contribution to final τ of all known systematics at 0.005, i.e., about 0.25σ. We conclude that we were unable to detect any systematic contribution to the 2015 Planckτ measurement as driven by low , and have limited it to well within our final statistical error budget.

Finally, since the submission of this paper, dedicated work on HFI data at low leads to a higher-precision determination of τ (Planck Collaboration Int. XLVI 2016) which is consistent with the one described in this paper. This latest work paves the way towards a future release of improved Planck likelihoods.

5.4.2. High- budget

We now turn to the high- likelihood. The approximate statistical model from which we build the likelihood function may turn out to be an unfaithful representation of the data for three main reasons. First the equations describing the likelihood or the parameters of those equations can be inaccurate. They are, of course, since we are relying on approximations, but we expect that in the regime where they are used our approximations are good enough not to bias the best fit or strongly alter the estimation of error bars. We call such errors due to a breakdown of the approximations a “methodological systematic”. We may also lump into this any coding errors. Second, our data model must include a faithful description of the relation between the sky and the data analysed, i.e., one needs to describe the transfer function and/or additive biases due to the non-ideal instrument and data processing. Again, an error in this model or in its parameters translates into possible errors that we call “instrumental systematics”. Finally, to recover the properties of the CMB, the contribution of astrophysical foregrounds must be correctly modelled and accounted for. Errors in this model or its parameters is denoted “astrophysical systematics”. When propagating each of these systematics to cosmological parameters, this is always within the framework of the ΛCDM model, as systematic effects can project differently into parameters depending on the details of the model.

We investigated the possibility of methodological systematics with massive Monte-Carlo simulations. One of the main technical difficulties of the high- likelihood is the computation of the covariance of the band powers. Appendix C.1.3 describes how we validated the covariance matrix, through the use of Monte-Carlo simulations, to better than a percent accuracy. This includes a first-order correction for the excess scatter due to point source masks, which can induce a systematic error in the covariance reaching a maximum of around 10% near the first peak and the largest scales (< 50), somewhat lower (about 5% or less) at other scales. In Sect. 3.6 we propagated the effect of those possible methodological systematics to the cosmological parameters and found a 0.1σ systematic shift on ns, when using the temperature data, which decreases when cutting the largest and most non-Gaussian modes. This is further demonstrated on the data in Sect. 5.1 where we vary the hybridization scale. At this stage it is unclear whether this is a sign of a breakdown of the Gaussian approximation at those scales, or if it is the result of the limitations of our point source correction to the covariance matrix. We did not try to correct for this bias in the likelihood and we assess this 0.1σ effect on ns to be the main contribution to the methodological systematics error budget.

Instrumental systematics are mainly assessed in three ways. First, given a foreground model, we estimate the consistency between frequencies and between the TT, EE and TE combinations at the spectrum and at the parameter level (removing some cross-spectra). For TT, the agreement is excellent, with shifts between parameters that are always compatible with the extra cosmic variance due to the removal of data when compared to the baseline solution (see Figs. 31 and 42). TE and EE inter-frequency tests reveal discrepancies between the different cross spectra that we assigned to leakage from temperature to polarization (see Fig. 40). In co-added spectra, these discrepancies tend to average out, leaving a few-μK2-level residual in the difference between the co-added TE and EE spectra and their theoretical predictions based on the TT parameters. Section 3.4.3 describes an effective model that succeeds in capturing some of that mismatch, in particular in TE. But as argued in Sect. 3.4.3 and Appendix C.3.5 one cannot, at this stage, use this model as-is to correct for the leakage, or to infer the level of systematic it may induce on cosmological parameters, due to a lack of a good prior on the leakage model parameters. However, cosmological parameters deduced from the current polarization likelihoods are in perfect agreement with those calculated from the temperature, within the uncertainty allowed by our covariance. The second way we assess possible instrumental systematics is by comparing the detset (DS) and the half-mission (HM) results. As argued in Sect. 3.4.4, the DS cross spectra are known to be affected by a systematic noise correlation that we correct for. Ignoring any uncertainty in this correction (which is difficult to assess), the overall shift between the HM- and DS-based parameters is of the order of 0.2σ (on ωb) at most on TT (Sect. 4.1.1 and Fig. 35), similar in TE and slightly worse in EE, particularily for ns. Since the uncertainty on the correlated noise correction is not propagated, those shifts are only upper bounds on possible instrumental systematics (at least those which would manifest differently in these two data cuts which are completely different as regards temporal systematics). Finally, in Sect. 3.7, we evaluate the propagation of all known instrumental effects to parameters. Due to the cost of the required massive end-to-end simulations, this test can only reveal large deviations; no such instrumental systematic bias is detected in this test. To summarize, our instrumental systematics budget is at most 0.2σ in temperature, slightly higher in EE, and there is no sign of bias due to temperature-to-polarization leakage that would not be compatible with our covariance (within the ΛCDM framework).

Finally, we assess the contribution of astrophysical systematics. Given the prior findings on polarisation, we only discuss the case of temperature here. The uncertainty on the faithfulness of the astrophysical model is relatively high, and we know from the DS/HM comparison that our astrophysical components certainly absorb part of the correlated noise that is not entirely captured by our model. In that sense, the recovered astrophysical parameters may be a biased estimate of the real astrophysical foreground contribution (due to the flexibility of the model which may absorb residual instrumental systematics provided they are sufficiently small). At small scales, the dominant astrophysical component is the point source Poisson term. We checked in Sect. 4.3 that the recovered point source contributions are in general agreement with models of their expected level. This is much less the case at 100 GHz and we argued in Sect. 4.3 that, nonetheless, an error in the description of the Poisson term is unlikely to translate into a bias in the cosmological parameters, as the point source contribution is negligible at all scales where the 100 GHz spectrum dominates the CMB solution. At large scales, the dust is our strongest foreground. We checked in Fig. 35 the effect of either marginalizing out the slope of the dust spectrum or removing the amplitude priors (i.e., making them arbitrarly wide). When marginalizing over the slope, one recovers a value compatible with the one in our model (−2.57 ± 0.038 whereas our model uses −2.63) and the cosmological parameters do not change (Sect. 4.1.2). When comparing the baseline likelihood result to CamSpec which uses a slightly different template we find a 0.16σ systematic shift in ns that can be attributed to the steeper dust template slope (−2.7) (Sect. 4.2). When ignoring the amplitude priors, a 0.2σ shift appears on ns (and As, due to its correlation with ns). However, in this case the level of dust contribution increases by about 20 μK2 in all spectra, which corresponds to more than doubling the 100 × 100 dust contribution. This level is completely ruled out by the 100 × 545 cross spectrum, which enables estimation of the dust contribution in the 100 GHz channel. The parameter shift can hence be attributed to a degeneracy between the dust model and the cosmological model broken by the prior on the amplitude parameters. We also use the fact that the dust distribution is anisotropic on the sky and evaluate the cosmological parameters on a smaller sky fraction. On TT there is no shift in the parameters that cannot be attributed to the greater cosmic variance on the smaller sky fraction. We are also making a simplifying assumption by describing the dust as a Gaussian field with a specific power spectrum. The numerical simulations (FFP9 and End-to-end) that include a realistic, anisotropic template for the dust contribution do not uncover any systematic effect due to that approximation. In the end, we believe that 0.2σ on ns is a conservative upper bound of our astrophysical systematic bias on the cosmological parameters. There is, however, a possibility of a residual instrumental bias affecting foreground parameters (but not cosmology), but we cannot, at this stage, provide quantitative estimates.

To summarize, our systematic error budget consists of a 0.1σ methodology bias on ns for TT, at most a 0.2σ instrumental bias on TT (on ωb), TE and possibly a slightly greater one on EE. The few-μK2-level leakage residual in polarization does not appear to project onto biases on the ΛCDM parameters. We conservatively evalute our astrophysical bias to be 0.2σ on ns. The astrophysical parameters might suffer from instrumental biases.

5.5. The low- “anomaly”

In Like13 we noted that the Planck 2013 low- temperature power spectrum exhibited a tension with the Planck best-fit model, which is mostly determined by high- information. In order to quantify such a tension, we performed a series of tests, concluding that the low- power anomaly was mainly driven by multipoles between = 20 and 30, which happen to be systematically low with respect to the model. The effect was shown to be also present (although less pronounced) using WMAP data (again, see Like13 and Page et al. 2007). The statistical significance of this anomaly was found to be around 99%, with slight variations depending on the Planck CMB solution or the estimator considered. This anomaly has drawn significant attention as a potential tracer of new physics (e.g., Kitazawa & Sagnotti 2015, 2014; Dudas et al. 2012; see also Destri et al. 2008), so it is worth checking its status in the 2015 analysis.

We present here updated results from a selection of the tests performed in 2013. While in Like13 we only concentrated on temperature, we now also consider low- polarization, which was not available as a Planck product in 2013. We first perform an analysis through the Hausman test (Polenta et al. 2005), modified as in Like13 for the statistic s1 = suprB(max,r), with max = 29 and H=CˆCVarCˆ,\begin{eqnarray} &&B(\ell_\mathrm{max},r) =\frac{1}{\sqrt{\ell_\mathrm{max}}}\sum_{\ell=2}^{\textrm{int}(\ell_\mathrm{max}r)}H_{\ell},\ \ \ r\in\left[0,1\right] \label{eq:hausman1} , \\ &&H_{\ell} =\frac{\hat{C_{\ell}}-C_{\ell}}{\sqrt{\textrm{Var}\,\hat{C_{\ell}}}} , \end{eqnarray}where Ĉ and C denote the observed and model power spectra, respectively. Intuitively, this statistic measures the relative bias between the observed spectrum and model, expressed in units of standard deviations, while taking the so-called “look-elsewhere effect” into account by maximizing s1 over multipole ranges. We use the same simulations as described in Sect. 2.3, which are based on FFP8, for the likelihood validation. We plot in Fig. 46 the empirical distribution for s1 in temperature and compare it to the value inferred from the PlanckCommander 2015 map described in Sect. 2 above. The significance for the Commander map has weakened from 0.7% in 2013 to 2.8% in 2015. This appears consistent with the changes between the 2013 and 2015 Commander power spectra shown in Fig. 2, where we can see that the estimates in the range 20 << 30 were generally shifted upwards (and closer to the Planck best-fit model) due to revised calibration and improved analysis on a larger portion of the sky. We also report in the lower panel of Fig. 46 the same test for the EE power spectrum, finding that the observed Planck low- polarization maps are anomalous only at the 7.7% level.

thumbnail Fig. 46

Top: empirical distribution for the Hausman s1 statistic for TT derived from simulations; the vertical bar is the observed value for the PlanckCommander map. Bottom: the empirical distribution of s1 for EE and the Planck 70 GHz polarization maps described in Sect. 2.

As a further test of the low- and high-Planck constraints, we compare the estimate of the primordial amplitude As and the optical depth τ, first separately for low and high multipoles, and then jointly. Results are displayed in Fig. 47, showing that the < 30 and the ≥ 30 data posteriors in the primordial amplitude are separated by 2.6σ, where the standard deviation is computed as the square root of the sum of the variances of each posterior. We note that a similar separation exists for τ, but it is only significant at the 1.5σ level. Fixing the value of the high- parameters to the Planck 2013 best-fit model slightly increases the significance of the power anomaly, but has virtually no effect on τ. A joint analysis using all multipoles retrieves best-fit values in As and τ which are between the low and high- posteriors. This behaviour is confirmed when the Planck 2015 lensing likelihood (Planck Collaboration XV 2016) is used in place of low- polarization.

thumbnail Fig. 47

Joint estimates of primordial amplitude As and τ for the data sets indicated in the legend. For low- estimates, all other parameters are fixed to the 2015 fiducial values, except for the dashed line, which uses the Planck 2013 fiducial. The PlanckTT+lowP estimates fall roughly half way between the low- and high- only ones.

Finally, we note a similar effect on Neff, which, in the high- analysis with a τ prior is about 1σ off the canonical value of 3.04, but is right on top of the canonical value once the lowP and its = 20 dip is included.

5.6. Compressed CMB-only high- likelihood

We extend the Gibbs sampling scheme described in Dunkley et al. (2013) and Calabrese et al. (2013) to construct a compressed temperature and polarization Planck high- CMB-only likelihood, Plik_lite , estimating CMB band-powers and the associated covariances after marginalizing over foreground contributions. Instead of using the full multi-frequency likelihood to directly estimate cosmological parameters and nuisance parameters describing other foregrounds, we take the intermediate step of using the full likelihood to extract CMB temperature and polarization power spectra, marginalizing over possible Galactic and extragalactic contamination. In the process, a new covariance matrix is generated for the marginalized spectra, which therefore includes foreground uncertainty. We refer to Appendix C.6.2 for a description of the methodology and to Fig. C.12 for a comparison between the multi-frequency data and the extracted CMB-only band-powers for TT, TE, and EE.

By marginalizing over nuisance parameters in the spectrum-estimation step, we decouple the primary CMB from non-CMB information. We use the extracted marginalized spectra and covariance matrix in a compressed, high-, CMB-only likelihood. No additional nuisance parameters, except the overall Planck calibration yP, are then needed when estimating cosmology, so the convergence of the MCMC chains is significantly faster. To test the performance of this compressed likelihood, we compare results using both the full multi-frequency likelihood and the CMB-only version, for the ΛCDM six-parameter model and for a set of six ΛCDM extensions.

We show in Appendix C.6.2 that the agreement between the results of the full likelihood and its compressed version is excellent, with consistency to better than 0.1σ for all parameters. We have therefore included this compressed likelihood, Plik_lite , in the Planck likelihood package that is available in the Planck Legacy Archive17.

5.7. Planck and other CMB experiments

5.7.1. WMAP-9

In Sect. 2.6 we presented the WMAP-9-based low- polarization likelihood, which uses the Planck 353 GHz map as a dust tracer, as well as the Planck and WMAP-9 combination. Results for these likelihoods are presented in Table 22, in conjunction with the Planck high- likelihood. Parameter results for the joint Planck and WMAP data set in the union mask are further discussed in Planck Collaboration XIII (2016) and Planck Collaboration XX (2016).

Table 22

Selected parameters estimated from Planck, WMAP, and their noise-weighted combination in low- polarization, assuming Planck in temperature at all multipoles.

thumbnail Fig. 48

Comparison of Planck and WMAP-9 CMB power spectra. Top: direct comparison. Noise spectra are derived from the half-ring difference maps. Bottom: residuals with respect to the PlanckΛCDM best-fit model. The error bars do not include the cosmic variance contribution (but the (brown) 1σ contour lines for the Likelihood best fit model do).

We now illustrate the state of agreement reached between the Planck 2015 data, in both the raw and likelihood processed form, and the final cosmological power spectra results from WMAP-9. In 2013 we noted that the difference between WMAP-9 and Planck data was mostly related to calibration, which is now resolved with the upward calibration shift in the Planck 2015 maps and spectra, as discussed in Planck Collaboration I (2016). This leads to the rather impressive agreement that has been reached between the two Planck instruments and WMAP-9.

Figure 48 (top panel) shows all the spectra after correction for the effects of sky masking, with different masks used in the three cases of the Planck frequency-map spectra, the spectrum computed from the Planck likelihood, and the WMAP-9 final spectrum. The Planck 70, 100, and 143 GHz spectra (which are shown as green, red, and blue points, respectively) were derived from the raw frequency maps (cross-spectra of the half-ring data splits for the signal, and spectra of the difference thereof for the noise estimates) on approximately 60% of the sky (with no apodization), where the sky cuts include the Galaxy mask, and a concatenation of the 70, 100, and 143 GHz point-source masks.

The spectrum computed from the Planck likelihood (shown in black as both individual and binned C values in Fig. 48) was described earlier in the paper. We recall that it was derived with no use of the 70 GHz data, but including the 217 GHz data. Importantly, since it illustrates the likelihood output, this spectrum has been corrected (in the spectral domain) for the residual effects of diffuse foreground emission, mostly in the low- range, and for the collective effects of several components of discrete foreground emission (including tSZ, point sources, CIB, etc.). This spectrum effectively carries the information that drives the likelihood solution of the Planck 2015 best-fit CMB anisotropy, shown in brown. Our aim here is to show the conformity between this Planck 2015 solution and the raw Planck data (especially at 70 GHz) and the WMAP-9 legacy spectrum.

The WMAP-9 spectrum (shown in magenta as both individual and binned C values) is the legacy product from the WMAP-9 mission, and it represents the final results of the WMAP team’s efforts to clean the residual effects of foreground emission from the cosmological anisotropy spectrum.

All these spectra are binned the same way, starting at = 30 with Δ = 40 bins, and the error-bars represent the error on the mean within each bin. In the low- range, especially near the first peak, the error calculation includes the cosmic variance contribution from the multipoles within each bin, which vastly exceeds any measurement errors (all the measurements shown here have high S/N over the first spectral peak), so we would expect good agreement between the errors derived for all the spectra in the completely signal-dominated range of the data.

The figure shows how WMAP-9 loses accuracy above ≈ 800 due to its inherent beam resolution and instrumental noise, and shows how the LFI 70 GHz data achieve improved fidelity over this range. HFI was designed to improve over both WMAP-9 and LFI in both noise performance and angular resolution, and the gains achieved are clearly visible, even over the relatively modest range of shown here, in the tiny spread of the individual C values of the Planck 2015 power spectrum. While the overall agreement of the various spectra, especially in the low- range, is noticeable in this coarse plot, it is also clear that the Planck raw frequency-map spectra do show excess power over the Planck best-fit spectrum at the higher end of the -range shown – the highest level at 70 GHz and the lowest at 143 GHz. This illustrates the effect of uncorrected discrete foreground residuals in the raw spectra.

A better view of these effects is seen in the bottom panel of Fig. 48. Here we plot the binned values from the top panel as deviations from the best-fit model. Naturally, the black bins of the likelihood output fit well, since they were derived jointly with the best-fit spectrum, while correcting for foreground residuals. The WMAP-9 points show good agreement, given their errors, with the Planck 2015 best fit, and illustrate very tight control of the large-scale residual foregrounds (at the low- range of the figure); beyond ~ 600 the WMAP-9 spectrum shows an increasing loss of fidelity. Planck raw 70, 100, and 143 GHz spectra show excess power in the lowest bin due to diffuse foreground residuals. The higher- range now shows more clearly the upward drift of power in the raw spectra, growing from 143 GHz to 70 GHz. This is consistent with the well-determined integrated discrete foreground contributions to those spectra. As previously shown in Planck Collaboration XXXI (2014, Fig. 8), the unresolved discrete foreground power (computed with the same sky masks as used here) can be represented in the bin near = 800 as levels of approximately 40 μK2 at 70 GHz, 15 μK2 at 100 GHz, and 5 μK2 at 143 GHz, in good agreement with the present figure.

5.7.2. ACT and SPT

Planck temperature observations are complemented at finer scales by measurements from the ground-based Atacama Cosmology Telescope (ACT) and South Pole Telescope (SPT). The ACT and SPT high-resolution data help Planck in separating the primordial cosmological signal from other Galactic and extragalactic emission, so as not to bias cosmological reconstructions in the damping-tail region of the spectrum. In 2013 we combined Planck with ACT (Das et al. 2014) and SPT (Reichardt et al. 2012) data in the multipole range 1000 << 10 000, defining a common foreground model and extracting cosmological parameters from all the data sets. Our updated “highL” temperature data include ACT power spectra at 148 and 218 GHz (Das et al. 2014) with a revised binning (Calabrese et al. 2013) and final beam estimates (Hasselfield et al. 2013), and SPT measurements in the range 2000 << 13 000 from the 2540 deg2 SPT-SZ survey at 95, 150, and 220 GHz (George et al. 2015). However, in this new analysis, given the increased constraining power of the Planck full-mission data, we do not use ACT and SPT as primary data sets. Using the same cuts as the 2013 analysis (i.e., ACT data at 1000 << 10 000 and SPT at > 2000) we only check for consistency and retain information on the nuisance foreground parameters that are not well constrained by Planck alone.

To assess the consistency between these data sets, we extend the Planck foreground model up to = 13 000 with additional nuisance parameters for ACT and SPT, as described in Planck Collaboration XIII (2016, Sect. 4). Fixing the cosmological parameters to the best-fit PlanckTT+lowP  base-ΛCDM model and varying the ACT and SPT foreground and calibration parameters, we find a reduced χ2 = 1.004 (PTE = 0.46), showing very good agreement between Planck and the highL data.

thumbnail Fig. 49

CMB-only power spectra measured by Planck (blue), ACT (orange), and SPT (green). The best-fit PlanckTT+lowP  ΛCDM model is shown by the grey solid line. ACT data at > 1000 and SPT data at > 2000 are marginalized CMB band-powers from multi-frequency spectra presented in Das et al. (2014) and George et al. (2015) as extracted in this work. Lower multipole ACT (500 << 1000) and SPT (650 << 3000) CMB power extracted by Calabrese et al. (2013) from multi-frequency spectra presented in Das et al. (2014) and Story et al. (2013) are also shown. The binned values in the range 3000 << 4000 appear higher than the unbinned best-fit line because of the binning (this is numerically confirmed by the residual plot in Planck Collaboration XIII 2016, Fig. 9).

As described in Planck Collaboration XIII (2016), we then take a further step and extend the Gibbs technique presented in Dunkley et al. (2013) and Calabrese et al. (2013; and applied to Planck alone in Sect. 5.6) to extract independent CMB-only band-powers from Planck, ACT, and SPT. The extracted CMB spectra are reported in Fig. 49. We also show ACT and SPT band-powers at lower multipoles as extracted by Calabrese et al. (2013). This figure shows the state of the art of current CMB observations, with Planck covering the low-to-high-multipole range and ACT and SPT extending into the damping region. We consider the CMB to be negligible at > 4000 and note that these ACT and SPT band-powers have an overall calibration uncertainty (2% for ACT and 1.2% for SPT).

The inclusion of ACT and SPT improves the full-mission Planck spectrum extraction presented in Sect. 5.6 only marginally. The main contribution of ACT and SPT is to constrain small components (e.g., the tSZ, kSZ, and tSZ×CIB) that are not well determined by Planck alone. However, those components are sub-dominant for Planck and are well described by the prior based on the 2013 Planck+highL solutions imposed in the Planck-alone analysis. The CIB amplitude estimate improves by 40% when including ACT and SPT, but the CIB power is also reasonably well constrained by Planck alone. The main Planck contaminants are the Poisson sources, which are treated as independent and do not benefit from ACT and SPT. As a result, the errors on the extracted Planck spectrum are only slightly reduced, with little additional cosmological information added by including ACT and SPT for the baseline ΛCDM model (see also Planck Collaboration XIII 2016, Sect. 4).

6. Conclusions

thumbnail Fig. 50

Planck 2015 CMB spectra, compared with the base ΛCDM fit to PlanckTT+lowP data (red line). The upper panels show the spectra and the lower panels the residuals. In all the panels, the horizontal scale changes from logarithmic to linear at the “hybridization” scale, = 29 (the division between the low- and high- likelihoods). For the residuals, the vertical axis scale changes as well, as shown by different left and right axes. We show \hbox{${\cal D}_\ell=\ell(\ell+1)C_\ell/(2\pi)$} for TT and TE, but C for EE, which also has different vertical scales at low- and high-.

The Planck 2015 angular power spectra of the cosmic microwave background derived in this paper are displayed in Fig. 50. These spectra in TT (top), TE (middle), and EE (bottom) are all quite consistent with the best-fit base-ΛCDM model obtained from TT data alone (red lines). The horizontal axis is logarithmic at < 30, where the spectra are shown for individual multipoles, and linear at ≥ 30, where the data are binned. The error bars correspond to the diagonal elements of the covariance matrix. The lower panels display the residuals, the data being presented with different vertical axes, a larger one at left for the low- part and a zoomed-in axis at right for the high- part.

The 2015 Planck likelihood presented in this work is based on more temperature data than in the 2013 release, and on new polarization data. It benefits from several improvements in the processing of the raw data, and in the modelling of astrophysical foregrounds and instrumental noise. Apart from a revision of the overall calibration of the maps, discussed in Planck Collaboration I (2016), the most significant improvements are in the likelihood procedures:

  • (i)

    a joint temperature-polarization pixel-based likelihood at ≤ 29, with more high-frequency information used for foreground removal, and smaller sky masks (Sects. 2.1 and 2.2);

  • (ii)

    an improved Gaussian likelihood at ≥ 30 that includes a different strategy for estimating power spectra from data-subset cross-correlations, using half-mission data instead of detector sets (which enables us to reduce the effect of correlated noise between detectors, see Sects. 3.2.1 and 3.4.3), and better foreground templates, especially for Galactic dust (Sect. 3.3.1) that lets us mask a smaller fraction of the sky (Sect. 3.2.2) and to retain large-angle temperature information from the 217 GHz map that was neglected in the 2013 release (Sect. 3.2.4).

We performed several consistency checks of the robustness of our likelihood-making process, by introducing more or less freedom and nuisance parameters in the modelling of foregrounds and instrumental noise, and by including different assumptions about the relative calibration uncertainties across frequency channels and about the beam window functions.

For temperature, the reconstructed CMB spectrum and error bars are remarkably insensitive to all these different assumptions. Our final high- temperature likelihood, referred to as “PlanckTT” marginalizes over 15 nuisance parameters (12 modelling the foregrounds, and 3 for calibration uncertainties). Additional nuisance parameters (in particular, those associated with beam uncertainties) were found to have a negligible impact, and can be kept fixed in the baseline likelihood. Detailed end-to-end simulations of the instrumental response to the sky analysed like the real data did not uncover hidden low-level residual systematics.

For polarization, the situation is different. Variation of the assumptions leads to scattered results, with greater deviations than would be expected due to changes in the data subsets used, and at a level that is significant compared to the statistical error bars. This suggests that further systematic effects need to be either modelled or removed. In particular, our attempt to model calibration errors and temperature-to-polarization leakage suggests that the TE and EE power spectra are affected by systematics at a level of roughly 1 μK2. Removal of polarization systematics at this level of precision requires further work, beyond the scope of this release. The 2015 high- polarized likelihoods, referred to as “Plik TE” and “Plik EE”, or “Plik TT, EE, TE” for the combined version, ignore these uncertain corrections. They only include 12 additional nuisance parameters accounting for polarized foregrounds. Although these likelihoods are distributed in the Planck Legacy Archive18, we stick to the PlanckTT+lowP choice in the baseline analysis of this paper and the companion papers such as Planck Collaboration XIII (2016), Planck Collaboration XIV (2016), and Planck Collaboration XX (2016).

We developed internally several likelihood codes, exploring not only different assumptions about foregrounds and instrumental noise, but also different algorithms for building an approximate Gaussian high- likelihood (Sect. 4.2). We compared these codes to check the robustness of the results, and decided to release:

  • (i)

    A baseline likelihood called Plik (available for TT, TE, EE, or combined observables), in which the data are binned in multipole space, with a bin-width increasing from Δ = 5 at ≈ 30 to Δ = 33 at ≈ 2500.

  • (ii)

    An unbinned version which, although slower, is preferable when investigating models with sharp features in the power spectra.

  • (iii)

    A simplified likelihood called Plik_lite in which the foreground templates and calibration errors are marginalized over, producing a marginalized spectrum and covariance matrix. This likelihood does not allow investigation of correlations between cosmological and foreground/instrumental parameters, but speeds up parameter extraction, having no nuisance parameters to marginalize over.

In this paper we have also presented an investigation of the measurement of cosmological parameters in the minimal six-parameter ΛCDM model and a few simple seven-parameter extensions, using both the new baseline Planck likelihood and several alternative likelihoods relying on different assumptions. The cosmological analysis of this paper does not replace the investigation of many extended cosmological models presented, e.g., in Planck Collaboration XIII (2016), Planck Collaboration XIV (2016), and Planck Collaboration XX (2016). However, the careful inspection of residuals presented here addresses two questions:

  • (i)

    a priori, is there any indication that an alternative model toΛCDM could provide a significantly better fit?

  • (ii)

    if there is such an indication, could it come from caveats in the likelihood-building (imperfect data reduction, foreground templates or noise modelling) instead of new cosmological ingredients?

Since this work is entirely focused on the power-spectrum likelihood, it can only address these questions at the level of 2-point statistics; for a discussion of higher-order statistics, see Planck Collaboration XVI (2016) and Planck Collaboration XVII (2016).

The most striking result of this work is the impressive consistency of different cosmological parameter extractions, performed with different versions of the Plik TT+tauprior or PlanckTT+lowP likelihoods, with several assumptions concerning: data processing (half-mission versus detector set correlations); sky masks and foreground templates; beam window functions; the use of two frequency channels instead of three; different cuts at low or high ; a different choice for the multipole value at which we switch from the pixel-based to the Gaussian likelihood; different codes and algorithms; the inclusion of external data sets like WMAP-9, ACT, or SPT; and the use of foreground-cleaned maps (instead of fitting the CMB+foreground map with a sum of different contributions). In all these cases, the best-fit parameter values drift by only a small amount, compatible with what one would expect on a statistical basis when some of the data are removed (with a few exceptions summarized below).

The cosmological results are stable when one uses the simplified Plik_lite likelihood. We checked this by comparing PlanckTT+lowP results from Plik and Plik_lite for ΛCDM, and for six examples of seven-parameter extended models.

Another striking result is that, despite evidence for small unsolved systematic effects in the high- polarization data, the cosmological parameters returned by the Plik TT, Plik TE, or Plik EE likelihoods (in combination with a τ prior or Planck lowP) are consistent with each other, and the residuals of the (frequency combined) TE and EE spectra after subtracting the temperature ΛCDM best-fit are consistent with zero. As has been emphasized in other Planck 2015 papers, this is a tremendous success for cosmology, and an additional proof of the predictive power of the standard cosmological model. It also suggests that the level of temperature-to-polarization leakage (and possibly other systematic effects) revealed by our consistency checks is low enough(on average over all frequencies) not to significantly bias parameter extraction, at least for the minimal cosmological model. We do not know yet whether this conclusion applies also to extended models, especially those in which the combination of temperature and polarization data has stronger constraining power than temperature data alone, e.g., dark matter annihilation (Planck Collaboration XIII 2016) or isocurvature modes (Planck Collaboration XX 2016). One should thus wait for a future Planck release before applying the Planck temperature-plus-polarization likelihood to such models. However, the fact that we observe a significant reduction in the error bars when including polarization data is very promising, since this reduction is expected to remain after the removal of systematic effects.

Careful inspection of residuals with respect to the best-fit ΛCDM model has revealed a list of anomalies in the Planck CMB power spectra, of which the most significant is still the low- temperature anomaly in the range 20 ≤ ≤ 30, already discussed at length in the 2013 release. In this 2015 release, with more data and with better calibration, foreground modelling, and sky masks, its significance has decreased from the 0.7% to the 2.8% level for the TT spectrum (Sect. 5.5). This probability is still small (although not very small), and the feature remains unexplained. We have also investigated the EE spectrum, where the anomaly, if any, is significant only at the 7.7% level.

Other “anomalies” revealed by inspection of residuals (and of their dependence on the assumptions underlying the likelihood) are much less significant. There are a few bins in which the power in the TT, TE, or EE spectrum lies 2–3σ away from the best-fit ΛCDM prediction, but this is not statistically unlikely and we find acceptable probability-to-exceed (PTE) levels. Nevertheless, in Sects. 3.8 and 4.1, we presented a careful investigation of these features, to see whether they could be caused by some imperfect modelling of the data. We noted that a deviation in the TT spectrum at ≈ 1450 is somewhat suspicious, since it is driven mostly by a single channel (217 GHz), and since it depends on the foreground-removal method. But this deviation is too small to be worrisome (1.8σ with the baseline Plik likelihood). As in the 2013 release, the data at intermediate would be fitted slightly better by a model with more lensing than in the best-fit ΛCDM model (to reduce the peak-to-trough contrast), but more lensing generically requires higher values of As and Ωch2 that are disfavoured by the rest of the data, in particular when high- information is included. This mild tension is illustrated by the preference for a value greater than unity for the unphysical parameter AL, a conclusion that is stable against variations in the assumptions underlying the likelihoods. However, AL is compatible with unity at the 1.8σ level when using the baseline PlanckTT likelihood with a conservative τ prior (to avoid the effect of the low- dip), so what we see here could be the result of statistical fluctuations.

This absence of large residuals in the Planck 2015 temperature and polarization spectra further establishes the robustness of the ΛCDM model, even with about twice as much data as in the Planck 2013 release. This conclusion is supported by several companion papers, in which many non-minimal cosmological models are investigated but no significant evidence for extra physical ingredients is found. The ability of the temperature results to pass several demanding consistency tests, and the evidence of excellent agreement down to the μK2 level between the temperature and polarization data, represent an important milestone set by the Planck satellite. The Planck 2015 likelihoods are the best illustration to date of the predictive power of the minimal cosmological model, and, at the same time, the best tool for constraining interesting, physically-motivated deviations from that model.


1

Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).

3

These rms estimates were computed with the PolSpice power-spectrum estimator (Chon et al. 2004) by averaging over 1000 noiseless simulations.

4

We assume here, and have checked in the data, that the noise-induced TQ and TU correlations are negligible.

5

The models considered have been derived by fixing all parameters except τ and As to their full multipole range 2015 best-fit values

7

To exactly mimic the procedure followed by the WMAP team, we exclude the signal correlation matrix from the noise component of the χ2 form. We have checked, however, that the impact of this choice is negligible for WMAP.

8

There is little point in comparing the scalings obtained for dust, as WMAP employs a model which is not calibrated to physical units.

10

We assume an -diagonal mixing matrix here. This is not necessarily the case, as sub-pixel beam effects, for example, can induce mode couplings. As discussed in Sect. 3.4.3, those were estimated in Like13 and found to be negligible for temperature. They are not investigated further in this paper.

11

Recall that these statements refer to the high- likelihood ( ≥ 30).

12

The definition of AL differs in PICO and CAMB ; see Appendix C.5.

13

The fitted values are 1.0000, 0.9999, 1.0000, 1.0000, 0.9987, 0.9986, 0.9992, 0.9989, 0.9989, 0.9981, 0.9989, 1.0000, and 0.9999 for detsets 100-ds1, 100-ds2, 143-ds1, 143-ds2, 143-5, 143-6, 143-7, 217-1, 217-2, 217-3, 217-4, 217-ds1, and 217-ds2, respectively.

14

During the revision of this paper, we noticed that the > 1000 case explores regions of parameter space that are outside the optimal PICO interpolation region, as also remarked by Addison et al. (2016). This inaccuracy mainly affected this particular test for constraints on ns and Ωbh2: the error bars for these parameters were underestimated by a factor of about 2 while the mean values were misestimated by about 0.8σ with respect to runs performed with CAMB . Nevertheless, we found that for all other parameters, and in all other likelihood tests presented in this section, this problem did not arise, since the explored parameter space was entirely contained in the PICO interpolation region so as to guarantee accurate results, as also detailed in Sect. C.5. Furthermore, this inaccuracy does not change any of the conclusions of this paper. We therefore decided to keep in Fig. 35 the results obtained with PICO   but we have added results for the > 1000 case obtained with CAMB (case “CAMB, lmin = 1000”).

15

These results were obtained with the PICO code, and are thus close to but not identical to those obtained with CAMB and reported in Planck Collaboration XIII (2016).

16

We note that the Hillipop mask was constructed partly so that Eq. (55) would be an accurate approximation. We find that for the radio contribution is is accurate to 2%, or 1σ, and for the dust contribution it is essentially exact.

19

We use the routine process_mask of the HEALPix package to obtain a map of the distance of each pixel of the mask from the closest null pixel. We then use a smoothed version of the distance map to build the Gaussian apodization. The smoothing of the distance map is needed to avoid sharp edges in the final mask.

20

Masking can in principle reduce the effective rank, but for the high sky fractions used in the Planck analysis, this is not an issue.

21

We remind the reader that, in this test, at each frequency we always use maxfreq=min(max,maxfreq,baseline)\hbox{$\lmax ^{\mathrm{freq}}=\mathrm{min}(\lmax ,\lmax ^{\mathrm{ freq,\,baseline} })$}, with maxfreq,baseline\hbox{$\lmax ^{\mathrm{freq,\, baseline}} $} the baseline max at each frequency as reported in Table 16 (e.g., in the max = 1404 case, we still use the 100 × 100 power spectrum through = 1197).

22

We discovered late in the preparation of this paper that in some of the tests the prior for the 143 × 217TE dust contamination was set inaccurately, with an offset of −0.3 μK2 at = 500. With our cuts, this spectrum contributes only at > 500 where the dust contamination is already small compared to the signal. We verified that this has no impact on the cosmology and on our conclusions.

23

Therefore, e.g., C143h2×217h1˜=(1+d143h2+d217h1)C143h2×217h1\hbox{$\tilde{C}_\ell^\mathrm{143h2~\times~217h1} = (1 + d_\mathrm{143h2} + d_\mathrm{217h1})\,C_\ell^\mathrm{143h2~\times~217h1}$}.

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, 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); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration.

We further 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 DPC infrastructure hosted at the Institut d’Astrophysique de Paris (France) and financially supported by CNES.

References

  1. Aad, G. E. A. 2014, Phys. Rev. D, 90, 052004 [NASA ADS] [CrossRef] [Google Scholar]
  2. Addison, G. E., Dunkley, J., & Spergel, D. N. 2012, MNRAS, 427, 1741 [NASA ADS] [CrossRef] [Google Scholar]
  3. Addison, G. E., Huang, Y., Watts, D. J., et al. 2016, ApJ, 818, 132 [NASA ADS] [CrossRef] [Google Scholar]
  4. Araujo, D., Bischoff, C., et al. QUIET Collaboration 2012, ApJ, 760, 145 [NASA ADS] [CrossRef] [Google Scholar]
  5. Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Audren, B., Lesgourgues, J., Benabed, K., & Prunet, S. 2013, J. Cosmol. Astropart. Phys., 2, 1 [Google Scholar]
  7. Battaglia, N., Natarajan, A., Trac, H., Cen, R., & Loeb, A. 2013, ApJ, 776, 83 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
  9. Benoît, A., Ade, P., Amblard, A., et al. 2003, A&A, 399, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJ, 757, L23 [NASA ADS] [CrossRef] [Google Scholar]
  11. Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 07, 034 [NASA ADS] [CrossRef] [Google Scholar]
  13. Bond, J. R., Contaldi, C., & Pogosyan, D. 2003, Roy. Soc. London Philosoph. Trans. Ser. A, 361, 2435 [Google Scholar]
  14. Calabrese, E., Hlozek, R. A., Battaglia, N., et al. 2013, Phys. Rev. D, 87, 103012 [NASA ADS] [CrossRef] [Google Scholar]
  15. Challinor, A., & Chon, G. 2005, MNRAS, 360, 509 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chon, G., Challinor, A., Prunet, S., Hivon, E., & Szapudi, I. 2004, MNRAS, 350, 914 [NASA ADS] [CrossRef] [Google Scholar]
  17. Chu, M., Eriksen, H. K., Knox, L., et al. 2005, Phys. Rev. D, 71, 103002 [NASA ADS] [CrossRef] [Google Scholar]
  18. Corasaniti, P. S., & Melchiorri, A. 2008, Phys. Rev. D, 77, 103507 [NASA ADS] [CrossRef] [Google Scholar]
  19. Crites, A. T., Henning, J. W., Ade, P. A. R., et al. 2015, ApJ, 805, 36 [NASA ADS] [CrossRef] [Google Scholar]
  20. Das, S., Marriage, T. A., Ade, P. A. R., et al. 2011, ApJ, 729, 62 [NASA ADS] [CrossRef] [Google Scholar]
  21. Das, S., Louis, T., Nolta, M. R., et al. 2014, JCAP, 04, 014 [Google Scholar]
  22. de Bernardis, P., Ade, P. A. R., Bock, J. J., et al. 2002, ApJ, 564, 559 [NASA ADS] [CrossRef] [Google Scholar]
  23. Delabrouille, J., Betoule, M., Melin, J.-B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Destri, C., de Vega, H. J., & Sanchez, N. G. 2008, Phys. Rev. D, 78, 023013 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dudas, E., Kitazawa, N., Patil, S. P., & Sagnotti, A. 2012, JCAP, 05, 012 [NASA ADS] [CrossRef] [Google Scholar]
  26. Dunkley, J., Komatsu, E., Nolta, M. R., et al. 2009, ApJS, 180, 306 [NASA ADS] [CrossRef] [Google Scholar]
  27. Dunkley, J., Calabrese, E., Sievers, J., et al. 2013, JCAP, 07, 025 [NASA ADS] [CrossRef] [Google Scholar]
  28. Durrer, R., Novosyadlyj, B., & Apunevych, S. 2003, ApJ, 583, 33 [NASA ADS] [CrossRef] [Google Scholar]
  29. Efstathiou, G. 2004, MNRAS, 349, 603 [NASA ADS] [CrossRef] [Google Scholar]
  30. Efstathiou, G. 2006, MNRAS, 370, 343 [NASA ADS] [CrossRef] [Google Scholar]
  31. Efstathiou, G., & Migliaccio, M. 2012, MNRAS, 423, 2492 [NASA ADS] [CrossRef] [Google Scholar]
  32. Elsner, F., & Wandelt, B. D. 2012, A&A, 542, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Eriksen, H. K., O’Dwyer, I. J., Jewell, J. B., et al. 2004, ApJS, 155, 227 [NASA ADS] [CrossRef] [Google Scholar]
  34. Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [NASA ADS] [CrossRef] [Google Scholar]
  35. Fendt, W. A., & Wandelt, B. D. 2007a, ArXiv e-prints [arXiv:0712.0194] [Google Scholar]
  36. Fendt, W. A., & Wandelt, B. D. 2007b, ApJ, 654, 2 [Google Scholar]
  37. Finelli, F., De Rosa, A., Gruppuso, A., & Paoletti, D. 2013, MNRAS, 431, 2961 [NASA ADS] [CrossRef] [Google Scholar]
  38. Fowler, J. W., Acquaviva, V., Ade, P. A. R., et al. 2010, ApJ, 722, 1148 [NASA ADS] [CrossRef] [Google Scholar]
  39. Galli, S., Benabed, K., Bouchet, F., et al. 2014, Phys. Rev. D, 90, 063504 [NASA ADS] [CrossRef] [Google Scholar]
  40. George, E. M., Reichardt, C. L., Aird, K. A., et al. 2015, ApJ, 799, 177 [NASA ADS] [CrossRef] [Google Scholar]
  41. Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
  42. Grainge, K., Carreira, P., Cleary, K., et al. 2003, MNRAS, 341, L23 [NASA ADS] [CrossRef] [Google Scholar]
  43. Gruppuso, A., de Rosa, A., Cabella, P., et al. 2009, MNRAS, 400, 463 [NASA ADS] [CrossRef] [Google Scholar]
  44. Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 103013 [NASA ADS] [CrossRef] [Google Scholar]
  45. Hanany, S., Ade, P. A. R., Balbi, A., et al. 2000, ApJ, 545, L5 [NASA ADS] [CrossRef] [Google Scholar]
  46. Hancock, S., & Rocha, G. 1997, in Microwave Background Anisotropies, eds. F. R. Bouchet, R. Gispert, B. Guiderdoni, & J. Trân Thanh Vân, 179–188 [Google Scholar]
  47. Hansen, F. K., Górski, K. M., & Hivon, E. 2002, MNRAS, 336, 1304 [NASA ADS] [CrossRef] [Google Scholar]
  48. Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1 [NASA ADS] [Google Scholar]
  49. Hasselfield, M., Moodley, K., Bond, J. R., et al. 2013, ApJS, 209, 17 [NASA ADS] [CrossRef] [Google Scholar]
  50. Hinshaw, G., Spergel, D. N., Verde, L., et al. 2003, ApJS, 148, 135 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hinshaw, G., Nolta, M. R., Bennett, C. L., et al. 2007, ApJS, 170, 288 [NASA ADS] [CrossRef] [Google Scholar]
  52. Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
  53. Hu, W., & White, M. 1997, Phys. Rev. D, 56, 596 [NASA ADS] [CrossRef] [Google Scholar]
  54. Hu, W., Sugiyama, N., & Silk, J. 1997, Nature, 386, 37 [NASA ADS] [CrossRef] [Google Scholar]
  55. Hu, W., Hedman, M. M., & Zaldarriaga, M. 2003, Phys. Rev. D, 67, 043004 [NASA ADS] [CrossRef] [Google Scholar]
  56. James, F., & Roos, M. 1975, Computer Phys. Commun., 10, 343 [NASA ADS] [CrossRef] [Google Scholar]
  57. Jones, W. C., Ade, P. A. R., Bock, J. J., et al. 2006, ApJ, 647, 823 [NASA ADS] [CrossRef] [Google Scholar]
  58. Jones, W. C., Montroy, T. E., Crill, B. P., et al. 2007, A&A, 470, 771 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Kamionkowski, M., Kosowsky, A., & Stebbins, A. 1997, Phys. Rev. D, 55, 7368 [NASA ADS] [CrossRef] [Google Scholar]
  60. Keisler, R., Reichardt, C. L., Aird, K. A., et al. 2011, ApJ, 743, 28 [NASA ADS] [CrossRef] [Google Scholar]
  61. Kilbinger, M., Benabed, K., Cappe, O., et al. 2011, ArXiv e-prints [arXiv:1101.0950] [Google Scholar]
  62. Kitazawa, N., & Sagnotti, A. 2014, JCAP, 04, 017 [NASA ADS] [CrossRef] [Google Scholar]
  63. Kitazawa, N., & Sagnotti, A. 2015, Modern Phys. Lett. A, 30, 1550137 [NASA ADS] [CrossRef] [Google Scholar]
  64. Knox, L., & Page, L. 2000, Phys. Rev. Lett., 85, 1366 [NASA ADS] [CrossRef] [Google Scholar]
  65. Kogut, A., Spergel, D. N., Barnes, C., et al. 2003, ApJS, 148, 161 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  66. Kovac, J. M., Leitch, E. M., Pryke, C., et al. 2002, Nature, 420, 772 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  67. Lesgourgues, J. 2011, ArXiv e-prints [arXiv:1104.2932] [Google Scholar]
  68. Lewis, A., & Bridle, S. 2002, Phys. Rev., D66, 103511 [Google Scholar]
  69. Mangilli, A., Plaszczynski, S., & Tristram, M. 2015, MNRAS, 453, 3174 [NASA ADS] [CrossRef] [Google Scholar]
  70. Mocanu, L. M., Crawford, T. M., Vieira, J. D., et al. 2013, ApJ, 779, 61 [NASA ADS] [CrossRef] [Google Scholar]
  71. Naess, S., Hasselfield, M., McMahon, J., et al. 2014, JCAP, 10, 07 [Google Scholar]
  72. Netterfield, C. B., Devlin, M. J., Jarosik, N., Page, L., & Wollack, E. J. 1997, ApJ, 474, 47 [NASA ADS] [CrossRef] [Google Scholar]
  73. Page, L., Nolta, M. R., Barnes, C., et al. 2003, ApJS, 148, 233 [NASA ADS] [CrossRef] [Google Scholar]
  74. Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
  75. Pearson, T. J., Mason, B. S., Readhead, A. C. S., et al. 2003, ApJ, 591, 556 [NASA ADS] [CrossRef] [Google Scholar]
  76. Peiris, H. V., Komatsu, E., Verde, L., et al. 2003, ApJS, 148, 213 [NASA ADS] [CrossRef] [Google Scholar]
  77. Percival, W. J., & Brown, M. L. 2006, MNRAS, 372, 1104 [NASA ADS] [CrossRef] [Google Scholar]
  78. Planck Collaboration. 2015, The Explanatory Supplement to the Planck 2015 results, http://wiki.cosmos.esa.int/planckpla/index.php/Main_Page (ESA) [Google Scholar]
  79. Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  82. Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Planck Collaboration XIII. 2014, A&A, 571, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  87. Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  88. Planck Collaboration XXXI. 2014, A&A, 571, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Planck Collaboration II. 2016, A&A, 594, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Planck Collaboration III. 2016, A&A, 594, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Planck Collaboration V. 2016, A&A, 594, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Planck Collaboration VI. 2016, A&A, 594, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  100. Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  101. Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  102. Planck Collaboration XIV. 2016, A&A, 594, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Planck Collaboration XV. 2016, A&A, 594, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  104. Planck Collaboration XVI. 2016, A&A, 594, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Planck Collaboration XVIII. 2016, A&A, 594, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Planck Collaboration XIX. 2016, A&A, 594, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Planck Collaboration XX. 2016, A&A, 594, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Planck Collaboration XXI. 2016, A&A, 594, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. Planck Collaboration XXIII. 2016, A&A, 594, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Planck Collaboration XXV. 2016, A&A, 594, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Planck Collaboration XXVIII. 2016, A&A, 594, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. Planck Collaboration Int. XVI. 2014, A&A, 566, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Planck Collaboration Int. XIX. 2015, A&A, 576, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  119. Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  120. Planck Collaboration Int. XLVI. 2016, A&A, submitted [Google Scholar]
  121. Polarbear Collaboration, Ade, P. A. R., Akiba, Y., et al. 2014, ApJ, 794, 171 [NASA ADS] [CrossRef] [Google Scholar]
  122. Polenta, G., Marinucci, D., Balbi, A., et al. 2005, JCAP, 11, 01 [Google Scholar]
  123. Pryke, C., Ade, P., Bock, J., et al. 2009, ApJ, 692, 1247 [NASA ADS] [CrossRef] [Google Scholar]
  124. Readhead, A. C. S., Myers, S. T., Pearson, T. J., et al. 2004, Science, 306, 836 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  125. Reichardt, C. L., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 694, 1200 [NASA ADS] [CrossRef] [Google Scholar]
  126. Reichardt, C. L., Shaw, L., Zahn, O., et al. 2012, ApJ, 755, 70 [NASA ADS] [CrossRef] [Google Scholar]
  127. Rocha, G., Contaldi, C. R., Bond, J. R., & Gorski, K. M. 2011, MNRAS, 414, 823 [NASA ADS] [CrossRef] [Google Scholar]
  128. Rocha, G., Contaldi, C. R., Colombo, L. P. L., et al. 2010, arXiv e-print [arXiv:1008.4948] [Google Scholar]
  129. Rosset, C., Tristram, M., Ponthieu, N., et al. 2010, A&A, 520, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Rudjord, Ø., Groeneboom, N. E., Eriksen, H. K., et al. 2009, ApJ, 692, 1669 [NASA ADS] [CrossRef] [Google Scholar]
  131. Scott, D., & White, M. 1994, in CMB Anisotropies Two Years after COBE: Observations, Theory and the Future, ed. L. M. Krauss, 214 [Google Scholar]
  132. Seljak, U. 1997, ApJ, 482, 6 [Google Scholar]
  133. Shaw, L. D., Rudd, D. H., & Nagai, D. 2012, ApJ, 756, 15 [NASA ADS] [CrossRef] [Google Scholar]
  134. Sievers, J. L., Achermann, C., Bond, J. R., et al. 2007, ApJ, 660, 976 [NASA ADS] [CrossRef] [Google Scholar]
  135. Souradeep, T., & Ratra, B. 2001, ApJ, 560, 28 [NASA ADS] [CrossRef] [Google Scholar]
  136. Spergel, D., Flauger, R., & Hlozek, R. 2015, Phys. Rev. D, 91, 023518 [NASA ADS] [CrossRef] [Google Scholar]
  137. Story, K. T., Reichardt, C. L., Hou, Z., et al. 2013, ApJ, 779, 86 [NASA ADS] [CrossRef] [Google Scholar]
  138. Taburet, N., Hernández-Monteagudo, C., Aghanim, N., Douspis, M., & Sunyaev, R. A. 2011, MNRAS, 418, 2207 [NASA ADS] [CrossRef] [Google Scholar]
  139. Tegmark, M. 1997, Phys. Rev. D, 55, 5895 [Google Scholar]
  140. Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001 [NASA ADS] [CrossRef] [Google Scholar]
  141. Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
  142. Trac, H., Bode, P., & Ostriker, J. P. 2011, ApJ, 727, 94 [NASA ADS] [CrossRef] [Google Scholar]
  143. Tristram, M., Macías-Pérez, J. F., Renault, C., & Santos, D. 2005a, MNRAS, 358, 833 [NASA ADS] [CrossRef] [Google Scholar]
  144. Tristram, M., Patanchon, G., Macías-Pérez, J. F., et al. 2005b, A&A, 436, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  145. Tucci, M., Toffolatti, L., de Zotti, G., & Martínez-González, E. 2011, A&A, 533, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  146. Wandelt, B. D., Hivon, E., & Górski, K. M. 2001, Phys. Rev. D, 64, 083003 [NASA ADS] [CrossRef] [Google Scholar]
  147. Wright, E. L., Bennett, C. L., Górski, K., Hinshaw, G., & Smoot, G. F. 1996, ApJ, 464, L21 [NASA ADS] [CrossRef] [Google Scholar]
  148. Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [Google Scholar]

Appendix A: Sky masks

This appendix provides details of the way we build sky masks for the high- likelihood. Since it is based on data at frequencies between 100 and 217 GHz, Galactic dust emission is the main diffuse foreground to minimize. We subtract the SMICA CMB temperature map (Planck Collaboration IX 2016) from the 353 GHz map and we adopt the resulting CMB-subtracted 353 GHz map as a tracer of dust. After smoothing the map with a 10deg Gaussian kernel, we threshold it to generate a sequence of masks with different sky coverage. Galactic masks obtained in this way are named B80 to B50, where the number gives the retained sky fraction fsky in percent (Fig. A.1).

For the likelihood analysis, we aim to find a trade-off between maximizing the sky coverage and having a simple, but reliable, foreground model of the data. The combination of masks and frequency channels retained is given in Table A.1. In order to get C-covariance matrices for the cosmological analysis that are accurate at the few percent level (cf. Sect. 3.5), we actually use apodized versions of the Galactic masks. The apodization corresponds to a Gaussian taper of width σ = 2deg19. Apodized Galactic masks are also used for the polarization analysis. The effective sky fraction of an apodized mask is fsky=iwi2Ωi/(4π)\hbox{$f_{\rm sky}=\sum_i w^2_i \Omega_i/(4\pi)$}, where wi is the value of the mask in pixel i and Ωi is the solid angle of the pixel.

All the HFI frequency channels, except 143 GHz, are also contaminated by CO emission from rotational transition lines. Here we are concerned with emission around 100 and 217 GHz, associated with the CO J = 1 → 0 and J = 2 → 1 lines, respectively. Most of the emission is concentrated near the Galactic plane and is therefore masked out by the Galactic dust masks. However, there are some emission regions at intermediate and low latitudes that are outside the quite small B80 mask we use at 100 GHz. We therefore create a mask specifically targeted at eliminating CO emission. The Type 3 CO map, part of the Planck 2013 product delivery (Planck Collaboration XIII 2014), is sensitive to low-intensity diffuse CO emission over the whole sky. It is a multi-line map, derived using prior information on line ratios and a multi-frequency component separation method. Of the three types of Planck CO maps, this has the highest S/N. We smooth this map with a σ = 120′ Gaussian and mask the sky wherever the CO line brightness exceeds 1KRJkms-1. The mask is shown in Fig. A.2, before apodization with a Gaussian taper of FWHM = 30′.

Finally, we include extragalactic objects in our temperature masks, both point sources and nearby extended galaxies. The nearby galaxies that are masked are listed in Table A.2, together with the corresponding cut radii. For point sources, we build conservative masks for 100, 143, and 217 GHz separately. At each frequency, we mask sources that are detected above S/N = 5 in the 2015 point-source catalogue (Planck Collaboration XXVI 2016) with holes of radius three times the σ=FWHM/ln8\hbox{$\sigma=\mathit{FWHM}{/}\!\sqrt{\ln8}$} of the effective Gaussian beam at that frequency. We take the FWHM values from the elliptic Gaussian fits to the effective beams (Planck Collaboration XXVI 2016), i.e., FWHM values of \hbox{$9\parcm66$}, \hbox{$7\parcm22$}, and \hbox{$4\parcm90$} at 100, 143, and 217 GHz, respectively. We apodize these masks with a Gaussian taper of FWHM = 30′. As already noted, these masks are designed to reduce the contribution of diffuse and discrete Galactic and extragalactic foreground emission in the “raw” (half-mission and detset) frequency maps used for the baseline high- likelihood.

thumbnail Fig. A.1

Unapodized Galactic masks B50, B60, B70, and B80, from orange to dark blue.

thumbnail Fig. A.2

Unapodized CO mask (fsky = 87%).

Table A.1

Galactic masks used for the high- analysis.

Table A.2

Masked nearby galaxies and corresponding cut radii.

The masks described in this appendix are used in the papers on cosmological parameters (Planck Collaboration XIII 2016), inflation (Planck Collaboration XX 2016), dark energy (Planck Collaboration XIV 2016), and primordial magnetic fields (Planck Collaboration XIX 2016), which are notable examples of the application of the high- likelihood. However, the masks differ from those adopted in some of the other Planck papers. For example, reconstructions of gravitational lensing (Planck Collaboration XV 2016) and integrated Sachs-Wolfe effect (Planck Collaboration XXI 2016), constraints on isotropy and statistics (Planck Collaboration XVI 2016), and searches for primordial non-Gaussianity (Planck Collaboration XVII 2016) mainly rely on the high-resolution foreground-reduced CMB maps presented in Planck Collaboration IX (2016). Those maps have been derived by four component-separation methods that combine data from different frequency channels to extract “cleaned” CMB maps. For each method, the corresponding confidence masks, for both temperature and polarization, remove regions of the sky where the CMB solution is not trusted. This is described in detail in Appendices AD of Planck Collaboration IX (2016). The masks recommended for the analysis of foreground-reduced CMB maps are constructed as the unions of the confidence masks of all the four component separation methods. Their sky coverages are fsky = 0.776 in temperature and fsky = 0.774 in polarization. Since component separation mitigates the foreground contamination even at relatively low Galactic latitudes, those masks feature a thinner cut along the Galactic plane than the ones described in this appendix. Nevertheless, propagation of noise, beam, and extragalactic foreground uncertainties in foreground-cleaned CMB maps is more difficult, and this is the main reason why we do not employ them in the baseline high- likelihood. We also note that the recommended mask for temperature foreground-reduced maps has a greater number of compact object holes than the masks used here. This is due to the fact that some component separation techniques can introduce contamination of sources from a wider range of frequencies than the approach considered here for the high- power spectra. According to the tests provided in Sect. C.1.4, such masks would result in sub-optimal performance of the analytic C-covariance matrices.

Appendix B: Low- likelihood supplement

Appendix B.1: Sherman-Morrison-Woodbury formula

In the Planck 2015 release we follow a pixel-based approach to the joint low- likelihood (up to = 29) of T, Q, and U. This approach treats temperature and polarization maps consistently at HEALPix resolution Nside = 16, as opposed to the WMAP low- likelihood, which incorporates polarization information from lower-resolution maps to save computational time (Page et al. 2007). The disadvantage of a consistent-resolution, brute-force approach lies in its computational cost (Like13), which may require massively parallel coding (and adequate hardware) in order to be competitive in execution time with the high- part of the CMB likelihood (see, e.g., Finelli et al. 2013 for one such implementation). Such a choice, however, would hamper the ease of code distribution across a community not necessarily specialized in massively parallel computing. Luckily, the Sherman-Morrison-Woodbury formula and the related matrix determinant lemma provide a means to achieve good timing without resorting to supercomputers. To see how this works, rewrite the covariance matrix from Eq. (3) in a form that explicitly separates the C to be varied from those that stay fixed at the reference model: M=XY=2cutCXYPXY+XY=cut+1maxCXY,refPXY+NXY=2cutCXYPXY+M0,\appendix \setcounter{section}{2} \begin{eqnarray} \tens{M} &=& \sum_{XY}\sum_{\ell =2}^{\ell_{\rm cut}} C_\ell^{XY} \tens{P}_\ell^{XY} +\sum_{XY} \sum_{\ell=\ell_{\rm cut}+1}^{\ell_{\rm max}} C_\ell^{XY,\rm ref} \tens{P}_\ell^{XY} +\tens{N}\\ &&\equiv \sum_{XY} \sum_{\ell =2}^{\ell_{\rm cut}} C_\ell^{XY} \tens{P}_\ell^{XY} +\tens{M}_0, \end{eqnarray}where we have effectively redefined the fixed multipoles as “high- correlated noise”, as far as the varying low- multipoles are concerned. Next, note that for fixed , PTT\hbox{$\tens{P}_{\ell}^{TT}$} has rank20λ = 2 + 1, and this matrix may therefore be decomposed as PTT=(VTT)TATTVTT\hbox{$\tens{P}_{\ell}^{TT} = (\tens{V}_\ell^{TT})^\tens{T}\, \tens{A}_\ell^{TT} \, \tens{V}_\ell^{TT}$}, where ATT\hbox{$\tens{A}_\ell^{TT}$} and VTT\hbox{$\tens{V}_\ell^{TT}$} are (λ × λ) and (λ × Npix) matrices, respectively, which depend only upon the unmasked pixel locations. A similar decomposition holds for the PEE,BB\hbox{$\tens{P}_{\ell}^{EE,BB}$} matrices, while PTE\hbox{$\tens{P}_{\ell}^{\textit{TE}}$} can be expanded in the [VTT,VEE]\hbox{$[\tens{V}_\ell^{TT},\tens{V}_\ell^{EE}]$} basis for the corresponding . We can then write M=VTA(C)V+M0,\appendix \setcounter{section}{2} \begin{equation} \tens{M} = \tens{V}^\tens{T} \tens{A}(C_\ell) \tens{V} +\tens{M}_0, \end{equation}(B.3)where V=[V2TT,V2EE,V2BB,...VcutBB]\hbox{$\tens{V} = [\tens{V}_2^{TT},\tens{V}_2^{EE},\tens{V}_2^{BB},\dots\tens{V}_{\ell_{\rm cut}}^{BB}]$} is an (nλ × Npix) matrix with nλ = 3 [(cut + 1)2−4], and A(C) is an (nλ × nλ) block-diagonal matrix (accounting for four modes removed in monopole and dipole subtraction). Each -block in the latter matrix reads [CTTATTCTEATE0CTEATECEEAEE000CBBABB].\appendix \setcounter{section}{2} \begin{equation} \begin{bmatrix} C_\ell^{TT} \tens{A}_\ell^{TT} & C_\ell^{\textit{TE}} \tens{A}_\ell^{\textit{TE}} & \tens{0} \\ C_\ell^{\textit{TE}} \tens{A}_\ell^{\textit{TE}} & C_\ell^{EE} \tens{A}_\ell^{EE} & \tens{0} \\ \tens{0} & \tens{0} & C_\ell^{BB} \tens{A}_\ell^{BB} \end{bmatrix}. \end{equation}(B.4)Finally, using the Sherman-Morrison-Woodbury identity and the matrix determinant lemma, we can rewrite the inverse and determinant of M as M-1=M0-1M0-1VT(A-1+VM0-1VT)-1VM0-1|M|=|M0||A||A-1+VM0-1VT|.\appendix \setcounter{section}{2} \begin{eqnarray} && \tens{M}^{-1} = \, \tens{M}_0^{-1} - \tens{M}_0^{-1} \tens{V}^\tens{T}(\tens{A}^{-1} + \tens{V} \tens{M}_0^{-1} \tens{V}^\tens{T})^{-1}\tens{V} {\tens{M}_0^{-1}} \\ && |\tens{M}| = \, |\tens{M}_0| \, |\tens{A}| \, |\tens{A}^{-1} + \tens{V} \tens{M}_0^{-1} \tens{V}^\tens{T}|~. \end{eqnarray}Because neither V nor M0 depends on C, all terms involving only their inverses, determinants, and products may be precomputed and stored. Evaluating the likelihood for a new set of C then requires only the inverse and determinant of an (nλ × nλ) matrix, not an (Npix × Npix) matrix. For the current data selection, described in Sects. 2.2 and 2.3, we find nλ = 2688, which is to be compared to Npix = 6307, resulting in an order-of-magnitude speed-up compared to the brute-force computation.

Appendix B.2: Lollipop

We performed a complementary analysis of low- polarization using the HFI data, in order to check the consistency with the LFI-based baseline result. The level of systematic residuals in the HFI maps at low is quite small, but comparable to the HFI noise (see Planck Collaboration VIII 2016), so these residuals should be either corrected, which is the goal of a future release, or accounted for by a complete analysis including parameters for all relevant systematic effects, which we cannot yet perform. Instead, we use Lollipop , a low- polarized likelihood function based on cross-power spectra. The idea behind this approach is that the systematics are considerably reduced in cross-correlation compared to auto-correlation.

At low multipoles and for incomplete sky coverage, the C statistic is not simply distributed and is correlated between modes. Lollipop uses the approximation presented in Hamimeche & Lewis (2008), modified as described in Mangilli et al. (2015) to apply to cross-power spectra. We restrict ourselves to the one-field approximation to derive a likelihood function based only on the EE power spectrum at very low multipoles. The likelihood function of the C given the data ˜C\hbox{$\tilde{C}_\ell$} is then 2lnP(C|˜C)=[Xg]T[Mf-1][Xg],\appendix \setcounter{section}{2} \begin{equation} -2\ln P(C_\ell|\tilde{C}_\ell)=\sum_{\ell \ell'} [X_g]^\tens{T}_\ell [M_f^{-1}]_{\ell \ell'} [X_g]_{\ell'}, \end{equation}(B.7)with the variable [Xg]=Cf+Og(˜C+OC+O)Cfid+O,\appendix \setcounter{section}{2} \begin{equation} \left[X_g\right]_\ell = \sqrt{ C_\ell^{f} + O_\ell} \,\, g{\left(\frac{\tilde{C}_\ell + O_\ell}{C_\ell + O_\ell}\right)} \,\, \sqrt{ C_\ell^{\rm fid} + O_\ell} , \end{equation}(B.8)where g(x)=2(xlnx1)\hbox{$g(x)=\sqrt{2(x-\ln x -1)}$}, Cfid\hbox{$C_\ell^{\rm fid}$} is a fiducial model and O is the offset needed in the case of cross-spectra. This likelihood has been tested on Monte Carlo simulations including both realistic signal and noise. In order to extract cosmological information on τ from the EE spectrum alone, we restrict the analysis to the cross-correlation between the HFI 100 and 143 GHz maps, which exhibits the lowest variance.

At large angular scales, the HFI maps are contaminated by systematic residuals coming from temperature-to-polarization leakage (see Planck Collaboration VIII 2016). We used our best estimate of the Q and U maps at 100 and 143 GHz, which we correct for residual leakage coming from destriping uncertainties, calibration mismatch, and bandpass mismatch, using templates as described in Planck Collaboration VIII (2016). Even though the level of systematic effects is thereby significantly reduced, we still have residuals above the noise level in null tests at very low multipoles (\hbox{$\ell \leqslant 4$}). To mitigate the effect of this on the likelihood, we restrict the range of multipoles to = 5−20.

Cross-power spectra are computed on the cleanest 50% of the sky by using a pseudo-C estimate (Xpol , an extension to polarization of the code described in Tristram et al. 2005a). The mask corresponds to thresholding a map of the diffuse polarized Galactic dust at large scales. In addition, we also removed pixels where the intensity of diffuse Galactic dust and CO lines is strong. This ensures that bandpass leakage from dust and CO lines does not bias the polarization spectra (see Planck Collaboration VIII 2016).

We construct the C correlation matrix using simulations including CMB signal and realistic inhomogeneous and correlated noise. In order to take into account the residual systematics, we derive the noise level from the estimated BB auto-spectrum where we neglect any possible cosmological signal. This over-estimates the noise level and ensures conservative errors. However, this estimate assumes by construction a Gaussian noise contribution, which is not a full description of the residuals.

We then sample the reionization optical depth τ from the likelihood, with all other parameters fixed to the Planck 2015 best-fit values (Planck Collaboration XIII 2016). Without any other data, the degeneracy between As and τ is broken by fixing the amplitude of the first peak of the TT spectrum (directly related to Ase− 2τ) at = 200. The resulting distribution is plotted in Fig. B.1. The best fit is at τ=0.064-0.016+0.015,zre=8.7-1.6+1.4,\appendix \setcounter{section}{2} \begin{equation} \tau=0.064_{-0.016}^{+0.015},\qquad z_\mathrm{re} = 8.7_{-1.6}^{+1.4} , \end{equation}(B.9)in agreement with the current Planck low- baseline (see Table 2), even though this result only relies on the EE spectrum between = 5 and 20.

thumbnail Fig. B.1

Distribution of the reionization optical depth τ using the Lollipop likelihood, based on the cross-correlation of the 100 and 143 GHz channels.

Appendix C: High- baseline likelihood: Plik

In this appendix, we provide detailed information on the Plik baseline likelihood used at high . First we describe in Sect. C.1 the Plik covariance matrix, by providing the equations we have implemented, by giving results from some of the numerical tests we carried out, and by describing our procedure to deal with the excess variance (as compared to the prediction of our approximate analytical model) due to the point source mask. Section C.2 validates the overall Plik implementation with Monte Carlo simulations of the full mission. For reference, Sect. C.3 gives the results of a large body of validation and stability tests on the actual data, including polarization in particular. We also discuss the numerical agreement of the temperature- and polarization-based results on base-ΛCDM parameters. Section C.4 describes how we calculate co-added CMB spectra from foreground-cleaned frequency power spectra. Section C.5 compares Plik cosmological results obtained using the PICO or CAMB codes. Finally, Sect. C.6 details how we marginalize over nuisance parameters to provide a fast but accurate CMB-only likelihood.

Appendix C.1: Covariance matrix

Appendix C.1.1: Structure of the covariance matrix

Here we summarize the mathematical formalism implemented to calculate the pseudo-power spectrum covariance matrices for temperature and polarization.

In the following, the fiducial power spectra C are assumed to be the smooth theory spectra multiplied by beam (b) and pixel window function (p) for detectors i and j, Ci,j=bibjp2(CCMB+CFG(fi,fj)),\appendix \setcounter{section}{3} \begin{equation} C_{\ell}^{i, j} = b_{\ell}^{i} \, b_{\ell}^{j} \, p_{\ell}^2 \left( C_{\ell}^{\mathrm{CMB}} + C_{\ell}^{\mathrm{FG}}(f_{i}, f_{j}) \right) , \end{equation}(C.1)where the fk denote the frequency dependence of the foreground contribution.

We now present the equations used to compute all the unique covariance matrix polarization blocks that can be formed from temperature and E-mode polarization maps (Hansen et al. 2002; Hinshaw et al. 2003; Efstathiou 2004; Challinor & Chon 2005; Like13). They approximate the variance of the biased pseudo-power spectrum coefficients, before correcting for the effects of pixel window function, beam, and mask.

TTTT block:

Var ( TT i,j , TT p,q ) C TT i,p C TT i,p C TT j,q C TT j,q Ξ TT ∅∅ , ∅∅ [ ( i,p ) TT , ( j,q ) TT ]