Issue 
A&A
Volume 594, October 2016
Planck 2015 results



Article Number  A9  
Number of page(s)  42  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201525936  
Published online  20 September 2016 
Planck 2015 results
IX. Diffuse component separation: CMB maps
^{1} APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{2} African Institute for Mathematical Sciences, 6−8 Melrose Road, Muizenberg, Cape Town, South Africa
^{3} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{4} AixMarseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326, 13388 Marseille, France
^{5} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{6} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{7} Atacama Large Millimeter/submillimeter Array, ALMA Santiago Central Offices, Alonso de Cordova 3107, Vitacura, Casilla 763 0355, Santiago, Chile
^{8} CGEE, SCS Qd 9, Lote C, Torre C, 4°; andar, Ed. Parque Cidade Corporate, CEP 70308200, Brasília, DF, Brazil
^{9} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{10} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{11} CRANN, Trinity College, Dublin, Ireland
^{12} California Institute of Technology, Pasadena, California, USA
^{13} Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{14} Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001 Teruel, Spain
^{15} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{16} Consejo Superior de Investigaciones Científicas (CSIC), Madrid, Spain
^{17} DSM/Irfu/SPP, CEASaclay, 91191 GifsurYvette Cedex, France
^{18} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{19} Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland
^{20} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{21} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, Oviedo, Spain
^{22} Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada
^{23} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{24} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{25} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA
^{26} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{27} Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, Florida, USA
^{28} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 0065 Helsinki, Finland
^{29} Department of Physics, Princeton University, Princeton, New Jersey, USA
^{30} Department of Physics, University of California, Santa Barbara, California, USA
^{31} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{32} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{33} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{34} Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, 00133 Roma, Italy
^{35} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16, 20133 Milano, Italy
^{36} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, 34127 Trieste, Italy
^{37} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00133 Roma, Italy
^{38} Discovery Center, Niels Bohr Institute, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{39} Discovery Center, Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{40} European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago, Chile
^{41} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{42} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{43} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{44} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{45} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 0065 Helsinki, Finland
^{46} INAF−Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35131 Padova, Italy
^{47} INAF−Osservatorio Astronomico di Roma, via di Frascati 33, Monte Porzio Catone, Italy
^{48} INAF−Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 34127 Trieste, Italy
^{49} INAF/IASF Bologna, via Gobetti 101, 40127 Bologna, Italy
^{50} INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy
^{51} INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{52} INFN, Sezione di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{53} INFN, Sezione di Roma 1, Università di Roma Sapienza, P.le Aldo Moro 2, 00185 Roma, Italy
^{54} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica, 1, 00185 Roma, Italy
^{55} INFN/National Institute for Nuclear Physics, via Valerio 2, 34127 Trieste, Italy
^{56} IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Grenoble Alpes, IPAG, CNRS, IPAG, 38000 Grenoble, France
^{57} IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, 411 007 Pune, India
^{58} ImperialCollege London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{59} Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{60} Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, 38000 Grenoble, France
^{61} Institut Universitaire de France, 103 bd SaintMichel, 75005 Paris, France
^{62} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{63} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98bis boulevard Arago, 75014 Paris, France
^{64} Institut für Theoretische Teilchenphysik und Kosmologie, RWTH Aachen University, 52056 Aachen, Germany
^{65} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{66} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0371 Oslo, Norway
^{67} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
^{68} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{69} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{70} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{71} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{72} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{73} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{74} Kazan Federal University, 18 Kremlyovskaya St., 420008 Kazan, Russia
^{75} LAL, Université ParisSud, CNRS/IN2P3, 91898 Orsay, France
^{76} LERMA, CNRS, Observatoire de Paris, 61 avenue de l’Observatoire, 75000 Paris, France
^{77} Laboratoire AIM, IRFU/Service d’Astrophysique − CEA/DSM − CNRS − Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{78} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{79} Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{80} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{81} Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{82} Lebedev Physical Institute of the Russian Academy of Sciences, Astro Space Centre, 84/32 Profsoyuznaya st., 117997 Moscow, Russia
^{83} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{84} McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8, Canada
^{85} National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{86} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{87} Niels Bohr Institute, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{88} Niels Bohr Institute, Copenhagen University, Blegdamsvej 17, 1165 Copenhagen, Denmark
^{89} Nordita (Nordic Institute for Theoretical Physics), Roslagstullsbacken 23, 106 91 Stockholm, Sweden
^{90} Optical Science Laboratory, University College London, Gower Street, London, UK
^{91} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{92} SMARTEST Research Centre, Università degli Studi eCampus, via Isimbardi 10, 22060 Novedrate (CO), Italy
^{93} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{94} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{95} Sorbonne UniversitéUPMC, UMR 7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{96} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{97} Space Sciences Laboratory, University of California, Berkeley, California, USA
^{98} Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167 KarachaiCherkessian Republic, Russia
^{99} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{100} The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics,Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{101} Theory Division, PHTH, CERN, 1211 Geneva 23, Switzerland
^{102} UPMC Univ. Paris 06, UMR 7095, 98bis boulevard Arago, 75014 Paris, France
^{103} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{104} Universities Space Research Association, Stratospheric Observatory for InfraredAstronomy, MS 23211, Moffett Field, CA 94035, USA
^{105} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{106} University of Granada, Instituto Carlos I de Física Teórica y Computacional, Granada, Spain
^{107} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received: 20 February 2015
Accepted: 23 April 2016
We present foregroundreduced cosmic microwave background (CMB) maps derived from the full Planck data set in both temperature and polarization. Compared to the corresponding Planck 2013 temperature sky maps, the total data volume is larger by a factor of 3.2 for frequencies between 30 and 70 GHz, and by 1.9 for frequencies between 100 and 857 GHz. In addition, systematic errors in the forms of temperaturetopolarization leakage, analoguetodigital conversion uncertainties, and very long time constant errors have been dramatically reduced, to the extent that the cosmological polarization signal may now be robustly recovered on angular scales ℓ ≳ 40. On the very largest scales, instrumental systematic residuals are still nonnegligible compared to the expected cosmological signal, and modes with ℓ< 20 are accordingly suppressed in the current polarization maps by highpass filtering. As in 2013, four different CMB component separation algorithms are applied to these observations, providing a measure of stability with respect to algorithmic and modelling choices. The resulting polarization maps have rms instrumental noise ranging between 0.21 and 0.27μK averaged over 55′ pixels, and between 4.5 and 6.1μK averaged over pixels. The cosmological parameters derived from the analysis of temperature power spectra are in agreement at the 1σ level with the Planck 2015 likelihood. Unresolved mismatches between the noise properties of the data and simulations prevent a satisfactory description of the higherorder statistical properties of the polarization maps. Thus, the primary applications of these polarization maps are those that do not require massive simulations for accurate estimation of uncertainties, for instance estimation of crossspectra and crosscorrelations, or stacking analyses. However, the amplitude of primordial nonGaussianity is consistent with zero within 2σ for all local, equilateral, and orthogonal configurations of the bispectrum, including for polarization Emodes. Moreover, excellent agreement is found regarding the lensing Bmode power spectrum, both internally among the various component separation codes and with the bestfit Planck 2015 Λ cold dark matter model.
Key words: cosmology: observations / polarization / cosmic background radiation / diffuse radiation
© ESO, 2016
1. Introduction
This paper, one of a set associated with the 2015 release of data from the Planck^{1} satellite, presents maps of the cosmic microwave background (CMB) anisotropies derived from the full Planck data set comprising 50 months of observations from the Low Frequency Instrument (LFI) and 29 months from the High Frequency Instrument (HFI; Planck Collaboration I 2016). This analysis updates the temperatureonly analysis of the first 15.5 months of Planck observations discussed in Planck Collaboration XII (2014), and presents the first CMB polarization maps derived from Planck observations.
Much of the Planck analysis effort since the 2013 data release has revolved around understanding and reducing instrumental systematic uncertainties. As summarized in Planck Collaboration I (2016), this work has been highly successful, reducing the net power from systematic errors in the HFI CMB channels by almost two orders of magnitude on large angular scales (Planck Collaboration XI 2016). The main contributions to these improvements have come from better temperaturetopolarization leakage modelling, reduced analoguetodigital conversion (ADC) errors, and improved modelling of very long time constants (VLTCs; Planck Collaboration VI 2016; Planck Collaboration VIII 2016). With these improvements, the Planck observations are now sufficiently free of instrumental artefacts to allow a robust determination of the CMB polarization anisotropies on intermediate and small angular scales, covering multipoles ℓ ≳ 20. However, as described both in this paper, and in Planck Collaboration VIII (2016) and Planck Collaboration XI (2016), residual systematics are still not negligible compared to the CMB signal on the very largest scales (ℓ ≲ 20), and these modes are therefore removed from the current maps using a highpass filter. The rate of progress is still excellent, however, and updated allscale maps with low largescale systematics are expected to be released in the near future.
For temperature, the most significant improvement in the Planck 2015 analysis pipeline is calibration based on the orbital CMB dipole rather than the solar dipole (Planck Collaboration I2016). This change, combined with a better understanding of both the Planck beams and transfer functions (Planck Collaboration VI 2016; Planck Collaboration VIII 2016), has reduced the uncertainties in gain calibration, and the agreement between LFI, HFI, and the Wilkinson Microwave Anisotropy Probe (WMAP) has improved to the level of the uncertainties. The calibration of the 70, 100, and 143 GHz channels determined from the CMB temperature power spectrum around the first acoustic peak agrees to within a few tenths of a percent. A similar comparison between 70 GHz and WMAP Vband, and between 100 GHz and WMAP W band shows agreement at the 0.5% level. For further details see Sect. 5.4, Table 4, and Fig. 6 of Planck Collaboration I (2016).
The component separation efforts of the Planck 2015 release are summarized in three papers. The current paper is dedicated to CMB extraction, and presents the main Planck 2015 CMB maps in both temperature and polarization. Planck Collaboration X (2016) addresses astrophysical component separation as implemented by Commander (Eriksen et al. 2004, 2008), a Bayesian component separation algorithm, and presents a global model of the microwave temperature sky ranging from 408 MHz to 857 GHz, including detailed maps of synchrotron, freefree, spinning dust, thermal dust, and CO emission. In addition, a few minor components (line emission around 90 GHz and the thermal SunyaevZeldovich effect near the Coma and Virgo clusters) are included in the model, as are instrumental parameters in the form of calibration, bandpasses, monopoles, and dipoles. The corresponding polarization model includes only synchrotron and thermal dust emission. Planck Collaboration XXV (2016) presents a detailed analysis of the foregrounds below about 100 GHz in both temperature and polarization. For detailed descriptions of the various foreground components relevant for microwave component separation, see either Sect. 2 of Planck Collaboration XII (2014) or Sect. 4 of Planck Collaboration X (2016).
The foreground amplitude relative to CMB polarization is such that effective foreground suppression is required for almost any cosmological analysis, but the optimal approach depends sensitively on the topic in question. For instance, because the fluctuation power of diffuse polarized foregrounds decays as a powerlaw in multipole moment ℓ (Page et al. 2007; Gold et al. 2011; Planck Collaboration X 2016; Planck Collaboration Int. XXX 2016), it is of greater importance for highℓ CMB power spectrum and likelihood estimation to minimize noise sensitivity and to marginalize over unresolved point sources, than it is to model diffuse foregrounds with high accuracy. In this case, it is more convenient to parameterize the residual foregrounds in terms of power spectrum models, and to marginalize over these in terms of a few global parameters, than to marginalize over a large number of perpixel foreground parameters. The Planck 2015 CMB likelihood therefore employs crossspectra coupled to simple harmonic space foreground modelling (Planck Collaboration XI 2016), rather than the detailed foreground modelling described in this paper.
Similarly, because of the low – but nonnegligible – level of residual instrumental systematics on large angular scales in the Planck 2015 data, the lowℓ polarization likelihood also implements a special purpose cleaning algorithm, in terms of a simple template fit including only the cleanest 30, 70, and 353 GHz channels (Planck Collaboration VI 2016; Planck Collaboration VIII 2016; Planck Collaboration XI 2016). This approach is similar to that adopted for the 9yr WMAP likelihood (Bennett et al. 2013), and allows for easy propagation of uncertainties from correlated noise in terms of full pixelpixel noise covariance matrices.
Other applications, however, such as gravitational lensing and integrated SachsWolfe (ISW) reconstructions (Planck Collaboration XV 2016; Planck Collaboration XXI 2016), constraints on isotropy and statistics (Planck Collaboration XVI 2016), searches for primordial nonGaussianity (Planck Collaboration XVII 2016), and constraints on global geometry and topological defects (Planck Collaboration XVIII 2016), require actual CMB maps, and these are all based on the products described in this paper.
As in 2013, we apply four complementary CMB component separation algorithms to the Planck 2015 sky maps. In alphabetical order, these are: (1) Commander, a parametric pixelbased Bayesian CMB Gibbs sampler (Eriksen et al. 2004, 2008); (2) NILC, a needletbased internal linear combination method (Basak & Delabrouille 2012, 2013); (3) SEVEM, which implements linear template fitting based on internal templates in pixel space (Leach et al. 2008; FernándezCobos et al. 2012); and (4) SMICA, a semiblind spectralmatching algorithm fully defined in harmonic space (Cardoso et al. 2008; Planck Collaboration XII 2014). These codes were all applied to CMB temperature reconstruction in Planck Collaboration XII (2014), and have now been extended to polarization, as described in Appendices A–D. In addition, each algorithm is applied to several subsets of the full data set, including halfring, halfmission, and yearly data splits (Planck Collaboration I 2016). Comparing the resulting maps, both between algorithms and data splits, provides a good understanding of instrumental and algorithmic uncertainties.
The paper is organized as follows. In Sect. 2 we describe the Planck 2015 data selection and preprocessing. In Sect. 3 we briefly review the component separation methods, deferring mathematical details to Appendices A–D. In Sect. 4 we present the derived CMB maps. In Sect. 5 we quantify the residual foreground emission present in the maps by crosscorrelation with foreground templates. In Sect. 6 we present angular power spectra and corresponding cosmological parameters. In Sects. 7 and 8 we consider higherorder statistics and gravitational lensing. In Sect. 9 we summarize the main features and limitations of these maps, and provide recommendations on their application. We conclude in Sect. 10.
2. Data selection and preprocessing
In this paper we use the fullmission Planck data (Planck Collaboration VI 2016; Planck Collaboration VIII 2016) and accompanying simulations (Planck Collaboration XII 2016). The CMB temperature maps are derived using all nine frequency channels, from 30 to 857 GHz. The CMB polarization maps are derived from the seven frequency channels sensitive to polarization, from 30 to 353 GHz. In most cases, the component separation methods use frequency channel maps as input, with the exception of the temperature analysis performed by Commander, which uses maps from subsets of detectors in each frequency channel, as specified in Table 1 of Planck Collaboration X (2016). Commander additionally uses the WMAP and the 408 MHz sky maps (Bennett et al. 2013; Haslam et al. 1982) for lowresolution temperature reconstruction, feeding into the lowℓ temperature likelihood; see Appendix A for further details.
The primary foregroundreduced CMB maps are derived from fullmission maps, maximizing signaltonoise ratio and minimizing destriping errors. In addition, a number of data splits are analysed to enable internal consistency checks and to make estimates of the properties of the CMB maps. These data splits can also be used in analyses where more than one map is required as input. For the purposes of component separation, each subset of the data must contain the same combination of frequency channels as the fullmission data set. The following data splits have been analysed:

maps from the first and second half of each pointing period (“half rings”; HR1 and HR2);

maps from odd and even years, consisting of year 1+3 maps for the LFI channels plus year 1 maps for the HFI channels, and year 2+4 maps for the LFI channels plus year 2 maps for the HFI channels (“years”; YR1 and YR2);

maps from the first and second half of the mission, consisting of year 1+2 maps for the LFI channels plus halfmission 1 maps for the HFI channels, and year 3+4 maps for the LFI channels plus halfmission 2 for the HFI channels (“half mission”; HM1 and HM2).
2.1. Data
The frequency maps are described in detail in Planck Collaboration II (2016) for LFI and Planck Collaboration VII (2016) for HFI. The maps have several features that are relevant to component separation. We summarize them here, and the reader is recommended to consult the references for further details.

Monopole and dipole contributions from the CMB, CIB, and other astrophysical components are estimated and removed during mapmaking (Planck Collaboration VI 2016; Planck Collaboration VIII 2016). This has an effect on component separation, and each method treats monopoles and dipoles in a different way, as described in Appendices A–D.

The HFI maps are corrected for zodiacal light emission (ZLE) by subtracting a model of the emission at the ring level during mapmaking (Planck Collaboration VIII 2016). This differs from the treatment in the 2013 component separation (Planck Collaboration XII 2014), where the ZLE model was not subtracted.

Leakage from intensity to polarization due to bandpass mismatches between detectors is estimated from the measurements of the bandpasses before launch (Planck Collaboration IV 2016; Planck Collaboration VI 2016; Planck Collaboration VIII 2016). These estimates are subtracted from the polarized maps.

The HFI maps at 100, 143, and 217 GHz are renormalized in order to correct for farsidelobe effects in the calibration (Planck Collaboration VII 2016).

Missing pixels are filled with the average values of pixels in the surrounding area. This area is defined as being within a radius of 1° for LFI maps and HFI frequency maps, and within a radius of 1.̊5 for HFI detector subset maps.
2.2. Simulations
To validate our results, we analyse realistic simulations of the Planck data set called the full focal plane 8 (FFP8) simulations. They are based on detailed models of the instrument and sky, and are described in full in Planck Collaboration XII (2016). We summarize their contents here.
2.2.1. CMB
The CMB was simulated using an input Λ cold dark matter (ΛCDM) model based on the 2013 cosmological parameter results (Planck Collaboration XVI 2014). The fiducial simulation contains no primordial tensor modes or primordial nonGaussianity. However, four variants of the same CMB realization have been produced that include nonzero values of the tensortoscalar ratio and nonGaussianity of a local type.
2.2.2. Foregrounds
The Planck sky model (PSM) has been used to simulate the foreground components. The intensity part of the simulation includes all astrophysical components that were identified in the 2013 release. The diffuse components that are relevant at low frequencies consist of synchrotron, freefree, and spinning dust. At high frequencies, CO, thermal dust, and CIB are included. The foreground modelling has been improved in this version of the simulations. In particular, Planck 353 GHz data were used to improve the frequency scaling of the dust emission. In polarization, the diffuse foregrounds are synchrotron and thermal dust. The extragalactic emission from radio and infrared sources has been simulated in intensity and polarization, and the SZ effect from clusters of galaxies has been included in intensity.
2.2.3. Simulated observations
Timeordered data (TOD) for each detector were simulated using the satellite pointing, and the individual detector beams, bandpasses, noise properties, and data flags. The same mapmaking used for the data is used to generate maps from the simulated TOD. All of the maps from subsets of the data have also been generated from the simulations.
Two versions of the maps are available, with and without bandpass mismatch leakage. The latter is simulated using the average bandpass for all detectors in a frequency channel, eliminating the leakage effect. The version of the maps without bandpass leakage is considered in this paper.
In addition to the fiducial maps, a set of 10 000 Monte Carlo realizations of CMB and noise has been generated. These realizations are intended to be used to assess the uncertainties on the results.
2.2.4. Mismatch between simulations and data
In analysing the simulations, a number of deficiencies have become evident. First, the amplitude of the CMB component does not match that of the data. This is the (expected) consequence of the fact that the CMB model for the simulations was specified before the recalibration of the Planck data between the 2013 release and the present one was completed. This mismatch can be mitigated by increasing the amplitude of the CMB simulations by 1.3% when comparing them to the data. Second, the noise properties of the simulated maps do not precisely match those of the data. This does not significantly affect the analysis of the CMB temperature maps, since they are signaldominated. However, it does affect the analysis of CMB polarization maps because they are more noisy. The noise mismatch appears to be scaledependent, since the adjustment of the amplitude of the noise simulations to match the data depends on the resolution of the maps. This is explored in Sect. 7.
Any analysis that relies on simulations to estimate the uncertainties of a result from the CMB polarization will be limited by these mismatches. Despite this, many analyses are possible, including those using crossspectrum, crosscorrelation, or stacking techniques.
3. Component separation methods
The four methods used by Planck to separate the CMB from diffuse foreground emission were described in detail in Planck Collaboration XII (2014). They are representative of the main approaches to component separation developed in recent years. The methods can be divided into two types. The first type assumes only knowledge of the blackbody spectrum of the CMB, and the foregrounds are removed by combining the multifrequency data to minimize the variance of the CMB component. The second type constructs an explicit parameterized model of the CMB and foregrounds with an associated likelihood, and the CMB component is obtained by maximizing or sampling from the posterior distribution of the parameters. Either type of method may be implemented in the map domain or in the harmonic domain. We recall briefly their main features and comment on their application to polarization data. Descriptions of the changes in each algorithm since 2013 are given in the appendices.

Commander (Eriksen et al. 2006, 2008) is a Bayesian parametric method that works in the map domain. Both the CMB and foregrounds are modelled using a physical parameterization in terms of amplitudes and frequency spectra, so the method is well suited to perform astrophysical component separation in addition to CMB extraction (Planck Collaboration X 2016). The joint solution for all components is obtained by sampling from the posterior distribution of the parameters given the likelihood and a set of priors. To produce a highresolution CMB map, the separation is performed at multiple resolutions with different combinations of input channels. The final CMB map is obtained by combining these solutions in the spherical harmonic domain. This obviates the need for the Ruler step that was used in 2013 to extend the Commander solution to high resolution. A lowresolution version of the separation is used to construct the temperature power spectrum likelihood for large angular scales, as described in Planck Collaboration XI (2016). We note that Commander employs detector and detector set maps rather than full frequency maps, and excludes some specific detector maps judged to have significant systematic errors. Thus, the selection of data is not identical between Commander and the other three methods. For further details regarding data selection and processing, see Appendix A.

NILC (Delabrouille et al. 2009) is an implementation of internal linear combination (ILC) that works in the needlet (wavelet) domain. The input maps are decomposed into needlets at a number of different angular scales. The ILC solution for the CMB is produced by minimizing the variance at each scale. This has the advantage that the weights used to combine the data can vary with position on the sky and also with angular scale. The solutions are then combined to produce the final CMB map.

SEVEM (FernándezCobos et al. 2012) is an implementation of the templatecleaning approach to component separation that works in the map domain. Foreground templates are typically constructed by differencing pairs of maps from the low and highfrequency channels. The differencing is done in order to null the CMB contribution to the templates. These templates are then used to clean each CMBdominated frequency channel by finding a set of coefficients to minimize the variance of the map outside of a mask. Thus SEVEM produces multiple foregroundcleaned frequency channel maps. The final CMB map is produced by combining a number of the cleaned maps in harmonic space.

SMICA (Cardoso et al. 2008) is a nonparametric method that works in the spherical harmonic domain. Foregrounds are modelled as a small number of templates with arbitrary frequency spectra, arbitrary power spectra, and arbitrary correlation between the components. The solution is obtained by minimizing the mismatch of the model to the auto and crosspower spectra of the frequency channel maps. From the solution, a set of weights is derived to combine the frequency maps in the spherical harmonic domain to produce the final CMB map. Maps of the total foreground emission in each frequency channel can also be produced. In the analysis performed for the 2013 release (Planck Collaboration XII 2014), SMICA was the method that performed best on the simulated temperature data.
3.1. Extension to polarization
The methods described above were applied to Planck temperature data for the 2013 release (Planck Collaboration XII 2014), and they have been extended to operate on polarization data in the present work. A key distinction between the methods is the choice of operating domain. Two of the methods, Commander and SEVEM, operate in the map domain, so it is most natural for them to perform the polarized component separation on the Q and U maps. The other two methods, NILC and SMICA, operate in the harmonic or needlet domain. An intrinsic part of the transform of polarized maps to these domains is the decomposition of Q and U into E and Bmodes, which is accomplished by using spherical harmonic transforms on the full sky. Thus these two methods perform their separation directly on E and B.
3.2. Outputs
In addition to producing CMB maps, each method provides “confidence” masks to define the region of the sky in which the CMB solution is trusted in temperature and polarization. The procedure each method uses to create the mask is described in Appendices A–D. The confidence masks are used to define masks for further analysis of the data. Two of the pipelines, Commander and SMICA also produce foreground products, described in Planck Collaboration X (2016).
The first 1000 Monte Carlo realizations of CMB and noise have been propagated through the four pipelines. This has been done twice, once using the parameters derived from the data, and once using the parameters derived from the fiducial FFP8 maps, to provide a set of simulations to accompany each data set.
Fig. 1 Preferred masks for analysing componentseparated CMB maps in temperature (left) and polarization (right). 
Fig. 2 Componentseparated CMB temperature maps at full resolution, FWHM 5′, N_{side} = 2048. 
The methods produce maps at different resolutions, as described in Appendices A−D. The products have been brought to a standard resolution to compare them and for distribution. Standard resolution temperature maps have a Gaussian beam of 5′ FWHM and HEALPix^{2} (Górski et al. 2005) resolution N_{side} = 2048. Standard resolution polarization maps have a Gaussian beam of 10′ FWHM and HEALPix resolution N_{side} = 1024. If maps are produced at higher resolution then they are downgraded to these standard resolutions. The downgrading procedure for maps is to decompose them into spherical harmonics (T, or E and B, as appropriate) on the full sky at the input HEALPix resolution. The spherical harmonic coefficients, a_{ℓm}, are convolved to the new resolution using (1)where b_{ℓ} is the beam transfer function, p_{ℓ} is the HEALPix pixel window function, and the “in” and “out” superscripts denote the input and output resolutions. They are then synthesized into a map directly at the output HEALPix resolution. Masks are downgraded in a similar way. The binary mask at the starting resolution is first downgraded like a temperature map. The smooth downgraded mask is then thresholded by setting pixels where the value is less than 0.9 to zero and all others to unity, to make a binary mask. This has the effect of enlarging the mask to account for the smoothing of the signal. In addition to the standard resolution products, the maps, masks, and Monte Carlo realizations have been downgraded to lower resolutions for analyses that need them, using the above procedure.
Fig. 3 Differences between the componentseparated CMB temperature maps from the 2013 and the 2015 releases. The maps have been smoothed to FWHM 80′ and downgraded to N_{side} = 128. 
The polarization maps and Monte Carlo realizations have been decomposed into E and Bmode maps, and downgraded to lower resolutions too, for analyses that work on E and B directly. The decomposition is done on the fullsky maps using spherical harmonic transforms. The CMB maps from the realspace methods, Commander and SEVEM, are inpainted before doing the decomposition. Both the standard Commander and SEVEM CMB maps are inpainted inside their corresponding confidence masks using a constrained realization (see Appendix A for details). For both methods, only the CMB maps are inpainted, not the Monte Carlo realizations, since this would be too computationally expensive. The other two methods, NILC and SMICA, work on E and Bmodes, so it is possible to make E and B maps directly from their outputs in addition the the standard Q and U maps.
4. CMB maps
In this section we present and discuss the componentseparated CMB maps in temperature and polarization. The maps are available to download from the Planck Legacy Archive (PLA)^{3}. For temperature, we compare the maps to those extracted from the nominal mission data in 2013. We also compare the maps from different methods to assess their consistency, and we use the FFP8 simulations to assess the accuracy of the methods and the robustness of the solutions. Throughout the discussion, we make use of appropriate masks in order to highlight differences at intermediate and high Galactic latitudes, ignoring the plane of the Galaxy, where differences are much higher due to the complexity of the foreground signal and its dominance over the CMB.
4.1. Temperature maps
Temperature confidence masks produced by the methods have been used to make combined masks for further analysis of the data. The first mask is constructed as the union of the Commander, SEVEM, and SMICA confidence masks. The NILC mask is not included in the union because it removes a significantly smaller fraction of the sky. This union mask has f_{sky} = 77.6%. We refer to it as UT78 and adopt it as the preferred mask for analysing the temperature maps. It is shown on the left in Fig. 1. An extended version of the mask has been constructed by adding to the UT78 mask those pixels where the the standard deviation between the four temperature maps is greater than 10 μK. This mask has f_{sky} = 76.1%, and we refer to it as UTA76. A union mask has been created for the FFP8 simulations in the same way as for the data. It has f_{sky} = 73.5%, and we refer to it as FFP8UT74.
Masks and statistics of componentseparated CMB maps from data and FFP8 simulations.
The CMB temperature maps produced by the four methods are shown in Fig. 2. No obvious differences are seen in these maps, and they appear visually consistent outside the mask. An important assessment of the robustness and consistency of the CMB T component separation solutions is provided in Fig. 3, which shows the differences between the Planck 2013 and 2015 maps for each method. Several interesting features may be seen in these differences, most of which correspond directly to a better understanding of the systematic uncertainties in the new maps.
Starting with Commander, the most striking features are largescale swaths tracing the Planck scanning strategy with a peaktopeak amplitude of approximately 10 μK. This pattern is very similar to that originally pointed out by Larson et al. (2015), who found this by subtracting the 9yr WMAP ILC map (Bennett et al. 2013) from a templatecleaned version of the Planck 2013 100 GHz map. Similar patterns are also seen in the Commander residual maps shown in Fig. 2 of Planck Collaboration X (2016), corresponding to detector maps that are rejected from the new 2015 analysis. These structures are primarily due to two effects, namely destriping errors from bandpass mismatch between detectors and far sidelobe contamination (Planck Collaboration VIII 2016; Planck Collaboration X 2016). By rejecting particularly susceptible channels in the updated Commander analysis, these errors are greatly reduced in the new map.
Turning to the other three difference maps in Fig. 3, we see that the residuals are internally very similar, but quite different from the Commander residuals. Clear traces of the ZLE are apparent, which is explained by the fact that it was not subtracted in the mapmaking in 2013, but is subtracted in the updated processing (Planck Collaboration VI 2014; Planck Collaboration VIII 2016). Commander, on the other hand, is less sensitive to residual ZLE, for the following two reasons. First, in 2013 it used channels only up to 353 GHz, which are less affected by the ZLE than the higher frequencies. Second, by virtue of fitting independent thermal dust spectral parameters (index and dust temperature) per pixel, it can efficiently absorb the ZLE in the thermal dust component. However, some of the ZLE may still be observed in the Commander differences; remnants of the “red arcs” typically seen in the second and fourth quadrant of the sky are just visible in the difference of the Commander solutions, while being very evident in all other cases. In addition, a quadrupole component aligned with the CMB dipole (Galactic coordinates (l,b) = (264°,48°); Planck Collaboration I 2016) can be seen. This is in part associated with the relativistic Doppler quadrupole. In 2013, the HFI data processing did not subtract this contribution, whereas the frequencyindependent component of the quadrupole is removed from the 2015 maps. The difference map thus contains the frequencyindependent part of this quadrupole, with small additional contributions due to frequencydependent residuals arising due to the change of frequency weights between 2013 and 2015. In a reference frame aligned with the CMB dipole, this corresponds to a 4–5 μK amplitude for the a_{20} component.
The residuals seen in Fig. 3 are small compared to the typical CMB anisotropies, with features mostly constrained to below 5–10 μK, with a distinct largescale pattern. In particular, these small differences are completely negligible for power spectrum and cosmological parameter estimation. The only cosmological application for which some care is needed is the study of large scale isotropy (Planck Collaboration XXIII 2014). An updated isotropy analysis of the new sky maps is presented in Planck Collaboration XVI (2016); no significant differences are found compared to the 2013 results. However, for studies of the CMB quadrupole−octopole alignment or 2point correlation function, correcting for the residual Doppler quadrupole may be desirable.
The residual from the frequencydependent component of the Doppler quadrupole has been estimated using the model described in Planck Collaboration XII (2016). The amplitude of the residual is found by propagating simulated maps of the frequencydependent part of the kinematic quadrupole through each of the component separation pipelines. The resulting maps give an estimate of the residual effect in the CMB maps after the removal of the frequencyindependent component, and may be applied to them as a further correction for the kinematic quadrupole effect. They have been made publicly available through the PLA as part of the ancillary data accompanying the CMB maps.
Table 1 summarizes the main properties of the CMB maps derived both from the data and from the fiducial set of FFP8 simulations. We evaluate standard deviations in two cases, corresponding to high (FWHM 10′, N_{side} = 1024) and intermediate resolution (FWHM 160′, N_{side} = 64). The values in parentheses are standard deviations of halfmission halfdifference (HMHD) maps, and they give an estimate of the level of uncertainties due to instrumental noise and systematic effects. The same quantities are given for the FFP8 simulations. At this level, results show good consistency for both data and simulations. The SEVEM maps have the lowest standard deviation, as measured by the HMHD maps, at small and intermediate angular scales.
Fig. 4 Pairwise difference maps between CMB temperature maps. As in the previous Fig. 3, the maps have been smoothed to FWHM 80′ and downgraded to N_{side} = 128. 
Fig. 5 Difference between output and input CMB temperature maps from FFP8 simulations. Smoothing and downgrading is performed as in Fig. 4. 
Fig. 6 Componentseparated CMB Q maps at resolution FWHM 10′, N_{side} = 1024. 
Pairwise differences between all four maps are shown in Fig. 4, after smoothing to 80′ FWHM and downgrading to N_{side} = 128. As expected, differences are largest close the edge of the mask, where the absolute foreground level is the highest. Comparing these with the corresponding maps from 2013 shown in figure 6 of Planck Collaboration XII (2014), and noting that the new colour bar spans a range that is 4 times narrower than the previous one, we see that the internal agreement between the four methods is substantially better in the new maps, typically by about a factor of 2.
Figure 5 shows the differences between the FFP8 outputs and the input CMB map. The residuals are smallest for NILC and largest for SEVEM. However, we note that the foreground model adopted for the simulations was chosen to be more complex than the real sky, in order to explicitly probe modelling errors. In particular, the simulated thermal dust frequency spectrum exhibits a strong positive (and spatially dependent) curvature at low frequency, which is neither captured in the parametric models adopted by Commander, nor easily modelled by the spatial templates adopted by SEVEM. This additional complexity makes it hard to draw any strong conclusions about the performance on the real data, for which the thermal dust spectrum may be very well approximated by a simple onecomponent modified blackbody with a nearly spatially constant spectral index (Planck Collaboration X 2016).
4.2. Polarization maps
We now turn our attention to the foregroundreduced CMB polarization maps. As discussed extensively in Planck Collaboration I (2016), Planck Collaboration VI (2016), and Planck Collaboration VIII (2016), the residual systematics in the Planck polarization maps have been dramatically reduced compared to 2013, by as much as two orders of magnitude on large angular scales. Nevertheless, on angular scales greater than 10°, corresponding to ℓ ≲ 20, systematics are still nonnegligible compared to the expected cosmological signal. Different combinations of input frequency channels for the component separation have been explored in order to mitigate the polarization residuals. However, it was not possible, for this data release, to fully characterize the largescale residuals from the data or from simulations. Therefore the CMB polarization maps provided in the current release have been highpass filtered to remove the large angular scales. This has been implemented by applying a cosine filter to the E and B spherical harmonic coefficients of the maps. This filter is defined as (2)and we have used ℓ_{1} = 20 and ℓ_{2} = 40. The same filtering has been applied to the FFP8 fiducial maps and to the Monte Carlo simulations.
Fig. 7 Componentseparated CMB U maps at resolution FWHM 10′, N_{side} = 1024. 
Fig. 8 20° × 20° patch of the highpass filtered Commander CMB polarization map, centred on the north ecliptic pole, (l,b) = (96°,30°). Each map is pixelized with a HEALPix resolution of N_{side} = 1024, and has an angular resolution of 10′ FWHM. The top row shows Q and U maps derived from the fullmission data set, the middle row shows the corresponding E and B maps, and the bottom row shows the E and Bmodes of the halfring halfdifference (HRHD) map. Note the characteristic “+” and “×” patterns in the Q and U maps, and the clear asymmetry between E and B in the full data set. Also note that the HRHD (null) E map is consistent with both the full and HRHD B maps. 
Fig. 9 Pairwise differences between CMB Q maps, after smoothing to FWHM 80′ and downgrading to N_{side} = 128. 
Fig. 10 Pairwise differences between CMB U maps, after smoothing and downgrading are performed as in Fig. 9. 
Fig. 11 Difference between output and input CMB Q maps from FFP8 simulations. Smoothing and downgrading are performed as in Figs. 9 and 10. 
Fig. 12 Difference between output and input CMB U maps from FFP8 simulations. Smoothing and downgrading are performed as in Fig. 11. 
As we did for temperature, individual polarization confidence masks derived for each method have been used to make combined masks for further analysis. Our first polarization mask is simply the union of the Commander, SEVEM, and SMICA confidence masks, which has f_{sky} = 77.6% and we refer to it as UP78. However, the 1point statistics analysis summarized in Sect. 7.1, revealed significant point source contamination in the Commander, NILC, and SMICA maps using this mask. SEVEM was not affected by this problem because it applies an inpainting technique to remove the brightest point sources (see Appendix C for further details). For this reason, two extended versions of the mask were created, the first by excluding in addition the pixels where the standard deviation between the CMB maps, averaged in Q and U, exceeds 4 μK. This mask has f_{sky} = 76.7%, and we refer to it as UPA77. The second was made by additionally excluding from the union mask the polarized point sources detected in each frequency channel. It has f_{sky} = 77.4%, and we refer to it as UPB77. This mask is shown on the right side of Fig. 1, and we adopt this as the preferred polarization mask, since it is physically better motivated than UPA77. Also, although it keeps a larger fraction of the sky than UPA77, it is sufficient to alleviate the point source contamination.
Masks have been made in a similar way for the FFP8 simulations. The union mask has f_{sky} = 76.3%, and we refer to it as FFP8UP76. An extended mask that also excludes polarized point sources has f_{sky} = 75.7%, and we refer to it as FFP8UPA76. FFP8UPA76 is the preferred mask for FFP8 polarization analysis.
Figures 6 and 7 show the highpass filtered Q and U Stokes parameter CMB maps after applying the UPB77 mask. The maps are shown at full resolution, and are thus dominated by instrumental noise, except in the regions at the ecliptic poles where integration time is greatest. Visually, the methods operating in the harmonic (NILC and SMICA) and spatial (Commander and SEVEM) domains are more similar to each other than the other methods.
In order to have better visual insight at the map level, in Fig. 8 we show a 20° × 20° patch of the highpassfiltered Commander polarization maps centred on the north ecliptic pole. The top row shows the fullmission Q and U maps. Note the characteristic “+” pattern in Q and “×” pattern in U; this is the expected signal for a pure E mode signal. To make this point more explicit, the middle row shows the same map decomposed into E and B components. There is a clear asymmetry between them, with E having visibly more coherent power than B, again as expected for an Edominated signal. Finally, the third row shows the halfring halfdifference (HRHD) E and B maps, illustrating the noise level in the fullmission maps. Comparing the middle and bottom rows, there is clearly an Emode excess in the fullmission map, whereas the corresponding fullmission Bmode map is consistent with the HRHD Bmode map. In addition, the HRHD Emode map is also consistent with both the fullmission and HRHD Bmode maps, suggesting that all are consistent with instrumental noise.
Pairwise differences between the four polarization maps are shown in Figs. 9 and 10. The NILC and SMICA solutions appear closest to each other. The regions of the sky that are most affected by differences appear to be those with a higher noise level, as may be seen by comparing to Figs. 6 and 7. For completeness, Figs. 11 and 12 show the differences between the FFP8 outputs and the input CMB map. The differences show a pattern similar to that of the noise, though with higher amplitude with respect to pairwise differences of solutions from data, possibly reflecting again the enhanced complexity of the simulated sky with respect to real data.
The combination of highpass filtering and noise makes visual comparison of these maps difficult. The rms summary provided in Table 1 is more informative in this respect. Comparing the HMHD rms values listed in parentheses, we see that NILC and SMICA have the lowest effective polarization noise levels at high angular resolution, with rms values that are roughly 20% lower than those observed for Commander and SEVEM. One plausible explanation for this difference is the angular resolution adopted for the fitting process by the four methods; whereas Commander and SEVEM perform the polarization analysis at 10′ FWHM resolution, NILC and SMICA adopt a 5′ FWHM resolution. When comparing the rms values at an angular resolution of 10′, as presented in Table 1, the maps from the latter two methods are smoothed by postprocessing to a lower resolution, whereas the maps from the two former codes are not.
This effect is not relevant at intermediate angular scales, for instance at 160′ FWHM, as shown in the bottom section of Table 1. On these angular scales, we see that the situation among the codes is reversed, and Commander provides a 20% lower effective noise than the other three codes.
5. Correlation with external templates
Correlation of the CMB maps with templates of foreground emission provides a first diagnostic of residual contamination in the maps. We compute the correlation coefficient r between two maps x and y as (3)where the index i runs over the N_{pix} pixels observed with the common mask, ⟨ x ⟩ = ∑ _{i}x_{i}/N_{pix}, σ_{x} = [ ∑ _{i}(x_{i}− ⟨ x ⟩ )^{2}/ (N_{pix}−1)] ^{1 / 2}, and similarly for y. Both maps and templates are smoothed to FWHM 1° and downgraded to N_{side} = 256 before computing the correlation, and the analysis is performed separately on temperature and polarization maps.
The foreground templates considered for temperature are: the 408 MHz radio survey of Haslam et al. (1982); the velocityintegrated CO map of Dame et al. (2001); the fullsky Hα template of Finkbeiner (2003); and the Planck 857 GHz channel map. The uncertainty in the value of r due to chance correlations between foregrounds and the cleaned maps is estimated by computing the correlations between the templates and the 1000 simulations of CMB and noise provided by each method. For polarization, the only template we consider is one constructed from the WMAP 9yr maps. The template is made by smoothing the K and Ka band maps to FWHM 1° and differencing them to remove the CMB contribution. This produces a template containing the lowfrequency polarized foreground emission and noise. The correlation analysis is done twice, once with the original maps and templates, and once with a highpass filtered version. The resulting coefficient factors and 1σ uncertainties are shown in Table 2.
Correlation coefficients of CMB maps with foreground templates for temperature and polarization.
For temperature, the results for all methods are compatible with no correlation within 1σ. From this, we conclude that there are no significant temperature foreground residuals with the same morphology as the templates in the map, to a precision set by cosmic variance. For polarization, the analysis of unfiltered maps shows that SEVEM and SMICA are compatible with no correlation at the 1σ level, Commander has a moderate level of correlation at the 2σ level, and for NILC we find a correlation around 4.5σ. For highpass filtered maps (labelled by HPF), the level of correlation is reduced for Commander and NILC to about 1.5σ, while it increases for SEVEM and SMICA to about 4 and 2.5σ, respectively. These levels may be due to accidental correlations of the map with the template that are not taken into account when computing the uncertainty. Neither the noise in the template nor the systematics in the maps at large angular scales are modelled, the latter being important for the unfiltered case. In the filtered case, the signal is reduced relative to the noise in the template, which could exacerbate spurious correlations. From this analysis, it is difficult to draw strong conclusions about the residual contamination in the polarized CMB maps.
6. Power spectra and cosmological parameters
In this section we evaluate the foregroundcleaned maps in terms of CMB power spectra and cosmological parameters. We emphasize that the Planck 2015 parameter results are not based on the highresolution foregroundcleaned CMB maps presented in this paper, but are instead derived from the likelihood described in Planck Collaboration XIII (2016). That likelihood combines the lowresolution Commander temperature map derived from Planck, WMAP, and 408 MHz with a templatecleaned LFI 70 GHz polarization map in a pixelbased lowℓ likelihood, and adopts a crossspectrumbased estimator for the highℓ temperature and polarization likelihood. The highℓ likelihood, called Plik, relies on a careful choice of masks along with templates and modelling, all in the power spectrum domain, to reduce the contribution from diffuse Galactic and discrete Galactic and extragalactic foreground emission. The templates and source models are marginalized over when estimating cosmological parameters.
The parameter estimation methodology that we use here is primarily a tool to evaluate the quality of the highresolution CMB maps and to assess overall consistency with the Planck 2015 likelihood. We start with the foregroundcleaned CMB maps and masks described in Sect. 4. The maps have been cleaned from diffuse foregrounds and, to a varying extent, from extragalactic foregrounds. We construct simplified templates for the residual extragalactic foregrounds that are marginalized over when estimating cosmological parameters.
While the Planck 2015 likelihood takes into account calibration and beam uncertainties, we have not done so in this analysis. Two of the methods, Commander and SMICA, fit for the relative calibration between frequency channels, but the uncertainties from this process are not propagated into the maps. The other two methods, NILC and SEVEM, assume that the frequency channels are perfectly calibrated. None of the four methods propagate the uncertainties in the beam transfer functions into the CMB maps.
We are interested in assessing the relative quality of the CMB maps, for which it is more important to assess the spread of parameters between methods and as a function of angular scale, rather than to provide absolute numerical values. However, since the CMB maps we describe here are the basis for the analysis of the statistical isotropy of the CMB, primordial nonGaussianity, and gravitational lensing by largescale structure, it is both important and reassuring that the parameter values that we find are reasonably close to those of Planck Collaboration XIII (2016).
Fig. 13 Power spectra of the the foregroundcleaned CMB maps. Left: TT power spectra evaluated using the UT78 mask. Right: EE power spectra evaluated using the UP78 mask. The thick lines show the spectra of the halfmission halfsum maps containing signal and noise. The thin lines show the spectra of the halfmission halfdifference maps, which give an estimate of the noise and some of the residual systematic effects. The black line shows the Planck 2015 bestfit CMB spectrum for comparison. 
Fig. 14 CMB TT and EE power spectra for each of the four foregroundcleaned maps. The first panel shows the TT bandpowers after subtracting the bestfit model of residual extragalactic foregrounds. The black points show the frequencymap based Plik spectrum, and the grey line shows the bestfit ΛCDM model from the Planck 2015 likelihood. The second panel shows the residuals of the TT bandpowers after subtracting the bestfit ΛCDM model. The third and fourth panels show the same information for the EE spectra. The bands are the same for all methods, but the points are offset horizontally for clarity. The parameters of the foreground model are marginalized over when estimating parameters (see Fig. 15). 
Fig. 15 Comparison of cosmological parameters estimated from the TT and EE spectra computed from the foregroundcleaned CMB maps. The labels on the horizontal axis give the spectrum and the ℓ_{max} used to obtain each set of results. The black points show the parameters obtained from the frequency mapbased Plik likelihood using the same spectrum and ℓ_{max}. For comparison, we also show the values of the parameters obtained with the full Planck 2015 likelihood (including multipoles up to ℓ_{max} = 2500) as the horizontal line surrounded by a grey band giving the uncertainties. The foreground model used for the CMB maps is the methodtailored fullsky model from the FFP8 simulations in TT and a single D_{ℓ} ∝ ℓ^{2} template in EE. 
Figure 13 shows the TT and EE power spectra of the foregroundcleaned maps, masked with an apodized versions of the UT78 mask in temperature and the UP78 mask in polarization. The power spectra of the halfmission halfsum (HMHS) data are shown as thick lines; those of the halfmission halfdifference (HMHD) data are shown as thin lines. The HMHS spectra contain signal and noise, whereas the HMHD spectra contain only noise and some residual systematic effects. The variations in the temperature noise spectra at low to intermediate multipoles are caused by the component separation methods optimizing the tradeoff between foreground signal and noise as a function of scale. At high multipoles, the same spectra are smooth because the relative weighting of the frequency channels is set by the noise levels. The breaks in the Commander noise spectra are caused by the hybridization of maps at different resolutions to make the final map.
Figure 14 shows the TT and EE power spectra of the CMB signal in the maps estimated using XFaster (Rocha et al. 2010, 2011) from the HMHS and HMHD maps. The first panel shows the four TT power spectra derived from the CMB maps, the Plik spectrum derived from the frequency channel maps, and the bestfit ΛCDM power spectrum derived from the Planck 2015 likelihood using multipoles up to ℓ = 2500. The second panel shows the TT residuals after subtraction of the bestfit theoretical spectrum. The third and fourth panels show the same information for the EE spectra.
The component separation pipelines are principally designed to remove diffuse foreground emission from the CMB, so residuals from unresolved compact sources remain in the maps at small angular scales. For the TT spectrum, these residuals are modelled by using the FFP8 simulations to construct extragalactic source templates as a function of multipole for each method individually, and marginalizing over the corresponding amplitudes during parameter estimation. These templates, however, are only as good as the inputs to the simulations. If the model of the extragalactic components is not correct, this will translate into an error in the shape of the templates. For the EE spectrum, we do not have a detailed model for the residuals, so we assume that it has a shape D_{ℓ} ∝ ℓ^{2}, and marginalize over its amplitude during cosmological parameter estimation. Figure 14 shows the power spectra after subtracting the bestfit model for high−ℓ extragalactic source residuals for TT and EE.
There are differences between the power spectra determined from the CMB maps and the spectrum determined by Plik from the frequency maps. The differences are most noticeable around the first peak in the TT power spectrum, where the CMB map spectra have deficits of power that range from 50–100μK^{2} compared to the Plik spectrum. The spectrum of the SEVEM map is the closest to the Plik spectrum.
Some of the differences can be ascribed to sample variance from the different sky coverage and the different combination of input data. As noted above, the temperature spectra from the CMB maps are determined using the UT78 mask, which has f_{sky} = 69% after apodization, whereas the Plik spectrum is determined from a combination of frequencychannel spectra using more conservative masks with f_{sky} = 66% at 100 GHz, f_{sky} = 57% at 143 GHz, and f_{sky} = 47% at 217 GHz. A similar situation obtains in polarization, where the CMB map spectra are determined using the UP78 mask with f_{sky} = 76% after apodization, but the Plik spectrum is obtained from the frequency maps using masks with f_{sky} = 70% at 100 GHz, f_{sky} = 50% at 143 GHz, and f_{sky} = 41% at 217 GHz.
Not all of the observed differences in the spectra can be explained by sample variance. There are differences between the spectra of the maps from the four component separation methods. Sample variance from sky coverage does not contribute to them because the spectra have all been determined using the same mask. One possible source of the differences is the relative calibration of the frequency channels in the component separation. If there is a mismatch in the assumed relative calibration between the channels, or an error in calibration is introduced in the internal recalibration by a method that does it, then the CMB map will not have the expected properties. The relative calibration of each frequency map propagates through the pipelines in an ℓdependent way, so, in general, an error in relative calibration produces an error in the effective beam transfer function of the CMB map. There is some evidence to support this hypothesis. During the revision of this paper it was found that the recalibration applied by SMICA to the 44 GHz channel was erroneously large. If the 44 GHz channel is not recalibrated relative to the other channels, then the power spectrum of the resulting SMICA map is much closer to that of the SEVEM map.
Another potential source of differences between the methods is the modelling of the diffuse foregrounds. The models used by Commander and NILC allow the properties of the foregrounds to vary across the sky, whereas SEVEM and SMICA assume they are fixed. It has been assumed that the component separation methods remove the diffuse foregrounds perfectly, leaving only the (compact) extragalactic foregrounds to be modelled in the parameter analysis. If this is not the case, then the level of the residuals will vary between methods, and it must be taken into account in the subsequent analysis.
Cosmological parameters from the componentseparated maps in both temperature and polarization are determined using XFaster power spectra, coupled to CosmoMC (Lewis & Bridle 2002) using a correlated Gaussian likelihood. Specifically, we include multipoles between ℓ_{min} = 50 and ℓ_{max}, where ℓ_{max} = 1000, 1500, or 2000 for temperature, and 1000 or 1500 for polarization. We adopt a standard 6parameter ΛCDM model, and, since lowℓ data are not used in the likelihood, impose an informative Gaussian prior of τ = 0.070 ± 0.006. As mentioned above, we construct foreground templates in TT for each CMB map by propagating the simulated fullsky FFP8 foreground maps through the respective pipeline and estimating the resulting power spectra, normalized to some pivotal multipole.
The resulting cosmological parameters are summarized in Fig. 15 for both TT and EE. Starting with the temperature cases, we first observe good overall internal agreement between the four component separation methods, with almost all differences smaller than 1σ within each multipole band. Second, we also observe good agreement with the bestfit Planck 2015 ΛCDM model derived from the likelihood, since most of the differences are within 1σ. The notable exception is the power spectrum amplitude, A_{s}e^{− 2τ}, which is systematically low at about the 2σ level for ℓ_{max} = 1000 for all methods.
In more detail, however, there is some evidence of internal tensions at the 1–2σ level, most clearly seen in Ω_{c}h^{2} and A_{s}e^{− 2τ}. For these two parameters, there are almost 1σ shifts for NILC, SEVEM, and SMICA going from ℓ_{max} = 1000 to 1500, and again from ℓ_{max} = 1500 to 2000. Commander appears to be somewhat more robust with respect to multipole range than the other three methods. Although some variation is indeed expected by statistical variation alone, the combination of the shapes of the power spectrum differences seen in Fig. 14 and the systematic parameter trends suggest systematic uncertainties at the 1−2σ level due to various unmodelled residuals, as discussed above. For this reason, we do not recommend using the componentseparated maps for cosmological parameter estimation at this time; for that purpose the Planck likelihood method is preferred, since it shows much better stability with respect to multipole range (Planck Collaboration XI 2016).
For polarization, the results are more ambiguous, with fluctuations relative to the temperature prediction beyond 2σ. Note, however, that the cosmic variance contributions to the temperature and polarization parameters are essentially independent, and the two estimates are therefore only expected to agree statistically, not pointbypoint. Nevertheless, SEVEM in particular appears to show evidence of larger deviations than expected in polarization for several parameters, and the θ parameter for Commander in polarization shows some hints of tension with respect to the corresponding temperature estimate.
The differences in the parameters are driven by the differences in the power spectra. In the absence of a model of errors in the relative calibration and the diffuse foreground residuals, they are absorbed in part by the extragalactic foreground model or by the cosmological parameters. However, the remaining differences are at a low level, so overall we find good consistency between the four different component separation methods in both temperature and polarization. In the next data release, we will investigate the effect of the relative calibration of the frequency channels and the modelling of the residual diffuse emission, and we will improve the extragalactic foreground model.
7. Higherorder statistics
We now consider the higherorder statistics of the CMB maps in the form of 1point statistics (variance, skewness, and kurtosis), Npoint correlation functions, and primordial nonGaussianity (f_{NL}). We focus in particular on the polarization maps and their degree of consistency with the FFP8 simulations. The temperature results are described in Planck Collaboration XVI (2016), Planck Collaboration XVII (2016), and Planck Collaboration XVIII (2016).
7.1. 1point statistics
The variance, skewness, and kurtosis of the maps, and the preprocessing steps needed to compute them, are described in detail in Monteserín et al. (2008), Cruz et al. (2011), and Planck Collaboration XVI (2016). The procedure is to normalize the data, d_{p}, in pixel p by its expected rms, . The rms is calculated from 1000 FFP8 realizations (Planck Collaboration XI 2016), for both CMB anisotropies and instrumental noise. To the extent that represents an accurate description of the data variance, and both the sky signal and the instrumental noise are Gaussiandistributed quantities, will be Gaussian distributed with zero mean and unit variance. In temperature, the variance of the instrumental noise is subtracted in order to determine the variance of the CMB signal in the data. Once the variance of the CMB signal is estimated, it is used to extract the skewness and kurtosis from the normalized map. This procedure is well established, and provides a direct test for the presence of residual foregrounds and of CMB Gaussianity.
The temperature analysis reveals that the Planck 2015 CMB maps and the FFP8 fiducial CMB maps are fully compatible with the Monte Carlo simulations, therefore they can be used for further statistical analyses. For more details about the temperature results and the FFP8 validation analysis please refer to Appendix E.2 and Planck Collaboration XVI (2016).
For polarization, since the Q and U maps are not rotationally invariant, we consider the polarized intensity . P is not Gaussiandistributed with zero mean, and its skewness and kurtosis are nonvanishing; however, by comparing the data with the Monte Carlo ensemble, we can test for residual foregrounds and quantify the performance of the component separation methods. Noise in the P maps is complicated; rather than trying to remove the noise contribution (as is done in the case of temperature), we compare the P maps with the Monte Carlo simulations of CMB and noise together. Table 3 gives the resulting lowertail probabilities, that is, the percentage of simulations that show a lower variance, skewness, or kurtosis than the P maps. This is done at three different resolutions (N_{side} = 1024, 256, and 64) for both the data maps (using mask UPB77) and the fiducial FFP8 maps (using mask FFP8UPA76).
Percentage of simulations showing a lower variance, skewness, or kurtosis than the P maps, the “lowertail probability”.
At lower resolutions, Table 3 shows good agreement between the fiducial CMB maps and the Monte Carlo simulations for all methods. At high resolution, skewness, and kurtosis are mostly in agreement with the simulations, with NILC and SMICA showing slightly high values for the kurtosis (lower tail probabilities of 98.5% and 99.4%, respectively). However, there is excess variance in the maps at high resolution, 2−3.5σ away from the mean of the simulations, with Commander deviating the least.
Using the individual components of the FFP8 fiducial maps (CMB, noise, and foregrounds), we are able to quantify the contributions to the statistics of the highresolution maps from each component separately. Figure 16 compares the variance (left column), skewness (middle column), and kurtosis (right column) extracted from the FFP8 Monte Carlo simulation P maps (histograms) with those extracted from the FFP8 fiducial realization of: the sum of CMB and noise (brown); the sum of CMB, noise, and thermal dust (light teal); the sum of CMB, noise, and radio point sources (violet); and the sum of CMB, noise, and all foregrounds (pink). The last case defines the lower tail probabilities in Table 3. We have also investigated other components, but found that their contributions are negligible.
We see that the only case that is fully compatible with the Monte Carlo ensemble for all methods is that of CMB and noise. In particular, adding the thermal dust component increases the variance slightly outside the acceptable range for both NILC and SMICA and the same effect becomes even stronger when adding the rest of the foreground components for these methods. Commander is only slightly affected by foregrounds by this measure, and remains within the acceptable range even for the full foreground model, while SEVEM is an intermediate case. For the kurtosis, we find relatively high sensitivity to radio source residuals (the violet and pink lines coincide), but very low sensitivity to diffuse foregrounds. We conclude that the anomalous statistics seen in analysis of the highresolution componentseparated FFP8 CMB maps are due to the foreground components, which could plausibly be caused by the additional complexity of the FFP8 foreground model with respect to the real sky. The Monte Carlo simulations are compatible with the CMB and noise components of the fiducial map, as they should be by design.
Fig. 16 Polarized intensity variance (left column), skewness (middle column), and kurtosis (right column) evaluated from the FFP8 Monte Carlo simulations (histogram) and from components of the fiducial FFP8 map at N_{side} = 1024 outside the FFP8UPA76 mask. The variance distributions have been normalized to the mean value of the Monte Carlo distributions for visualization purposes. Coloured vertical lines correspond to different combinations of components: the sum of CMB and noise is shown in brown; the sum of CMB, noise, and thermal dust is shown in light teal; the sum of CMB, noise, and radio point sources is shown in violet; the sum of CMB, noise, and all foregrounds is shown in pink. 
Fig. 17 Polarized intensity variance evaluated from the FFP8 Monte Carlo simulations (histogram) and from the Planck 2015 maps (vertical pink lines) outside the UPB77 mask. Columns from left to right show different resolutions (N_{side} = 1024, 256, and 64), while rows show results for the four component separation methods. Unlike in Fig. 16, the variance distributions are not normalized here. 
Figure 17 shows the variance of the Monte Carlo simulations (from which the lower tail probabilities are summarized in the top half of Table 3) compared with that of the Planck 2015 CMB P maps. The maps analysed are the sum of CMB and noise signal. Since the CMB is practically the same for all methods, the different average values of the variance of the Monte Carlo simulations are a reflection of somewhat higher noise on small angular scales in the SEVEM and Commander maps than in the SMICA and NILC maps. Looking at the data (pink bars), we immediately see that they do not match the simulations. No simulation has a variance as high as the data for any method at any resolution. As already mentioned, this discrepancy is due to an underestimation of the noise in the FFP8 simulations (Planck Collaboration XII 2016). In the variance column of Table 3, in parentheses, we report the excess variance of the data in percent relative to the mean variance of the Monte Carlo simulations. The excess is 3−4% at N_{side} = 1024, increasing to 10−20% at N_{side} = 64. Moreover, we see that there are some extreme values for both the skewness and kurtosis, in particular at low resolution for Commander and at high resolution for NILC. However, as long as the variance in the simulations is inaccurate, it is difficult to decide whether residual foregrounds, a true nonGaussian feature in the map, or noise underestimation is the cause. For now, we simply conclude that the simulations are inconsistent with the data, and this effectively prevents studies of higherorder statistics of this kind for the polarization maps. For temperature, the same is not true because of the much higher signaltonoise ratio, although care is warranted even then when probing into the noisedominated regime above ℓ ≈ 1500−2000.
We also performed the 1point analysis on the Planck 2015 polarization maps using the UP78 mask. In this case, the values of the skewness and kurtosis for SMICA, NILC, and Commander were significantly affected by the presence of point sources. Conversely, the results for SEVEM were consistent with the expected distribution because this method applies an inpainting technique to remove the signal from the brightest point sources (see Appendix C for more details). This motivated the construction of the UPB77 mask, which excluded the brightest point sources detected in polarization, and strongly alleviated this problem for the other three methods, as can be seen in the results in Table 3. The same behaviour was also found for the FFP8 fiducial maps.
7.2. Npoint correlation functions
Realspace Npoint correlation functions are a useful diagnostic of the statistics of CMB maps, complementary to harmonic analyses. In this section we describe their application to the Planck 2015 CMB polarization maps. Results for the FFP8 CMB maps are given in Appendix E.3. Details of their application to temperature maps may be found in Planck Collaboration XVI (2016).
For observed fields X measured in a fixed relative orientation on the sky, the Npoint correlation function is defined as (4)where the unit vectors span an Npoint polygon on the sky. Assuming statistical isotropy, Npoint functions are functions only of the geometrical configuration of the Npoint polygon. In the case of the CMB, the fields X correspond to ΔT and two Stokes parameters Q and U describing the linear polarization of the radiation in direction . In standard CMB conventions, Q and U are defined with respect to the local meridian of the spherical coordinate system of choice. However, Q and U form a spin2 field and depend on a rotational coordinate system transformation (Leahy et al. 2010). To obtain coordinatesystemindependent Npoint correlation functions the Stokes parameters are rotated with respect to a local coordinate system defined by the centre of mass of the polygon (see Gjerløw et al. 2010 for details). The Stokes parameters in this new “radial” system are denoted by Q_{r} and U_{r}.
Given rotationally invariant quantities X ∈ {ΔT,Q_{r},U_{r}}, the correlation functions are estimated by simple product averages over all sets of N pixels fulfilling the geometric requirements set by θ_{1},...,θ_{2N−3}, which characterize the shape and size of the polygon, (5)The pixel weights are introduced to reduce noise or mask boundary effects. Masks set weights to 1 for included pixels and 0 for excluded pixels.
The shapes of the polygons selected for this analysis are pseudocollapsed and equilateral configurations for the 3point function, and a rhombic configuration for the 4point function, comprising two equilateral triangles sharing a common side. The 4point function is only computed in the analysis of temperature maps. We use the same definition of pseudocollapsed as in Eriksen et al. (2005), that is, an isosceles triangle where the length of the baseline falls within the second bin of the separation angles. The length of the longer edges of the triangle, θ, parameterizes its size. Analogously, in the case of the equilateral triangle and rhombus, the size of the polygon is parameterized by the length of the edge, θ. We note that these functions are chosen because of ease of implementation, not because they are better suited for testing Gaussianity than other configurations. In the following, all results refer to the connected 4point function.
We analyse the Planck 2015 CMB temperature and highpass filtered polarization maps at resolution FWHM 160′, N_{side} = 64. We used the downgraded version of the UP78 mask in the analysis.
The Npoint functions are used to test the quality of the CMB estimates derived from the Planck data, and are shown in Fig. 18. We show the differences between the Npoint functions for the highpass filtered CMB maps and the corresponding mean values estimated from the 1000 Monte Carlo simulations. The probabilities of obtaining values of the χ^{2} statistic for the Planck fiducial ΛCDM model at least as large as for the CMB map are given in Table 4. The results of the analysis of the temperature maps can be found in Planck Collaboration XVI (2016).
Fig. 18 Difference between the Npoint functions for the highpass filtered N_{side} = 64Planck 2015 CMB estimates and the corresponding means estimated from 1000 Monte Carlo simulations. The Stokes parameters Q_{r} and U_{r} were locally rotated so that the correlation functions are independent of coordinate frame. The first row shows results for the 2point function, from left to right, TQ_{r}, TU_{r}, Q_{r}Q_{r}, and U_{r}U_{r}. The second row shows results for the pseudocollapsed 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}, and the third row shows results for the equilateral 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r} and U_{r}U_{r}U_{r}. The red solid, orange tripledotdashed, green dashed and blue dotdashed lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using SMICA simulations. See Sect. 7.2 for the definition of the separation angle θ. 
The results for the data deviate significantly from the Monte Carlo simulations for almost all Npoint functions involving at least two polarization fields. The smallest deviation is seen for the Commander map. The Npoint functions for Monte Carlo simulations have smaller variance than for the data. By comparing HMHD maps with Monte Carlo simulations of noise corresponding to the CMB estimates, we find that the amplitude of the noise is underestimated by 18%. After adjusting the Monte Carlo simulations to compensate, the Npoint functions are more consistent with data. However, it is difficult to draw conclusions about the Gaussianity of the polarization maps from these results, since the Npoint functions themselves are used to estimate the mismatch between the noise Monte Carlo simulations and the data.
7.3. Primordial nonGaussianity
Primordial nonGaussianity is often measured in terms of the amplitude, , of the quadratic corrections to the gravitational potential, as well as by means of the 3point correlation function based on different triangle configurations. The results from these calculations for the foregroundcleaned CMB maps are presented in Planck Collaboration XVII (2016). Compared to the previous data release, we can now include both temperature and polarization bispectra in the analysis. We thus consider bispectrum f_{NL} estimates obtained from temperature and polarization data only, as well as the full constraint from all eight possible TTT, TTE, EET, and EEE combinations.
Results obtained from application of the KomatsuSpergelWandelt (KSW) estimator (Komatsu et al. 2005) to the CMB maps after subtraction of the lensingISW correlation (see below) are listed in Table 5 for various geometrical configurations that have been considered (Planck Collaboration XVII 2016). It is interesting to evaluate the impact of polarization data on the estimation of f_{NL} measurements. By considering only the temperature bispectrum for the SMICA map, we obtain , while the polarization alone yields . Uncertainties are evaluated by means of Gaussian FFP8 simulations. We find consistency between pipelines at the 1σ level. The results confirm and extend to polarization the absence of evidence of nonGaussianity of primordial origin estimated through the bispectrum. The performance of the methods is also tested using Gaussian and nonGaussian FFP8 simulations, showing that SMICA and SEVEM give the results closest to the inputs in both cases (see Sect. 7.3 of Planck Collaboration XVII 2016).
An interesting case to consider is that of the ISWlensing bispectrum (Planck Collaboration XXI 2016). Unlike with primordial shapes, there is a specific prediction for the expected amplitude of the ISWlensing threepoint signal in a given cosmological model. We can therefore use this shape to check whether the expected level of NG is recovered, verifying that different component separation methods neither spuriously add nor remove any signal, at least in the squeezed limit where this shape is peaked. This is indeed the case. By normalizing the ISWlensing shape in such a way as to have an f_{NL} amplitude of 1 in the bestfit model, we recover f_{NL}(SMICA) = 0.85 ± 0.2 from the full analysis including all bispectra, f_{NL}(SMICA) = 0.6 ± 0.3 from temperature alone, and f_{NL}(SMICA) = 4.7 ± 6.0 from polarization alone.
Amplitude of primordial nonGaussianity, f_{NL}, estimated by the KSW estimator.
The constraint from polarization alone is, as expected, much looser than the one from temperature alone. However, polarization does not have a negligible impact on the final combined measurement, mostly due to contributions coming from TTE configurations. Adding polarization reduces the final uncertainty by about 30%, as well as moving the recovered amplitude parameter closer to its expected value of 1. The implications of these results in terms of the physics of the early Universe, as well as the study of many additional shapes, are discussed in Planck Collaboration XVII (2016), which explains the algorithmic details and all results from the procedure summarized here, and in Planck Collaboration XX (2016), which discusses the implications for inflationary physics.
8. Gravitational lensing by largescale structure
Gravitational lensing from intervening matter imprints a nonGaussian signature in the CMB temperature and polarization maps, which in turn can be exploited to extract the gravitational potential integrated along the line of sight back to the surface of last scattering. Planck accurately measures the lensing potential over most of the sky (Planck Collaboration XV 2016). As it remaps the CMB polarization, lensing partially transforms primordial Emodes into Bmodes, resulting in a secondary Bmode spectrum peaking at around ℓ = 1000. By forming a weighted product of the Emodes from componentseparated CMB maps and the reconstructed lensing potential, it is possible to generate a map of the expected lensinginduced Bmodes of CMB polarization. Crosscorrelating this lensing Bmode template with the total observed Bmodes provides an indirect measurement of the lensing Bmode power spectrum (Planck Collaboration PIP116, in prep.; Planck Collaboration XV 2016). Lensing Bmode template maps were synthesized for the four component separation methods considered in this paper, using the Stokes parameter maps and the lensing potential reconstruction from the intensity map, as discussed in Planck Collaboration PIP116 (in prep.). A common mask was generated by combining the union polarization mask and the lensing potential 80% mask of Planck Collaboration PIP116 (in prep.). After apodization using a cosine function over 3°, the resulting mask preserves an effective sky fraction of about 60%.
Fig. 19 Lensinginduced Bmode power spectra in the componentseparated polarization CMB maps. The solid line represents the bestfit cosmology from the Planck data release in 2015. Error bars were evaluated using a semianalytical approximation validated over the FFP8 simulations, as described in Planck Collaboration PIP116 (in prep.). 
The lensing Bmode power spectra for Commander, NILC, SEVEM, and SMICA, which have been measured by crosscorrelating with the corresponding foreground cleaned polarization maps, are shown in Fig. 19. The four CMB polarization solutions lead to consistent lensing Bmode power spectrum measurements within 1σ over the entire probed multipole range. In addition, the χ^{2} relative to the 2015 PlanckΛCDM base model (Planck Collaboration XIII 2016) of the lensing Bmode band powers in eight multipole bins are 12.2, 9.0, 8.0, and 7.0, with corresponding probabilitytoexceed values of 14, 34, 44, and 54% for Commander, NILC, SEVEM, and SMICA, respectively, which indicate the agreement of the Planck lensing Bmode signal with theoretical expectations. We measure amplitude fits with respect to the Planck base model of
corresponding to 10, 11, 10, and 12σ detections of lensing Bmodes for Commander, NILC, SEVEM, and SMICA, respectively.
9. Summary and recommendations
We now summarize the results, and provide a critical analysis of the applicability of the derived maps for cosmological purposes.
Starting with the temperature case, we have shown that the four CMB maps are in excellent agreement overall. The amplitudes of the pairwise difference maps are smaller than 5 μK over most of the sky on large angular scales, and the highℓ power spectra agree to around 1σ. Correspondingly, the differences with respect to the Planck 2013 maps are typically smaller than 10 μK over most of the sky, and the morphology of the differences is well understood in terms of improved treatment of systematic errors in the 2015 analysis. We conclude that the Planck 2015 temperature maps provide a more accurate picture of the CMB sky than the 2013 temperature maps, having both lower noise and lower levels of systematics, and we expect these maps to find the same cosmological applications as the previous generation of maps. However, we emphasize that these maps are not cleaned of highℓ foregrounds, such as extragalactic point source or SZ emission. Residuals lead to biases in cosmological parameters at the 2σ level, beyond ℓ ≃ 2000. Cosmological analyses using small angular scales must therefore take care to marginalize over such foregrounds, as appropriate.
For polarization, the situation is significantly more complicated, due to two different problems. First, a low level of residual systematics on large angular scales prevents a faithful CMB polarization reconstruction for multipoles ℓ ≲ 20. These modes are therefore removed by highpass filtering in the current maps. Any cosmological analysis of these maps must take into account the corresponding transfer function in order to avoid biases. Second, due to the current noise mismatch between the FFP8 simulations and the data, we strongly caution against using the polarization maps provided here for any cosmological analysis that depends sensitively on the assumed noise level. Nevertheless, the maps should prove useful for a number of other applications that do not require detailed noise simulations, for instance estimation of crossspectra, crosscorrelations, or stacking analyses, and we therefore release the maps to the public despite these limitations.
Considering the four component separation methods in greater detail, the results may be distinguished according to two criteria, namely data selection and basis functions. While Commander performs data selection at the detector (or detectorset) map level, rejecting potentially problematic maps, the other three methods employ frequency channel maps, and thereby maximize the signaltonoise ratio. Likewise, while Commander and SEVEM perform their analyses in pixel space, NILC and SMICA perform all operations in harmonic space. These distinctions can explain many of the qualitative differences discussed in the previous sections.
We make the following recommendations regarding the use of the four maps. First and foremost, we strongly recommend that any cosmological analysis based on these maps consider all four maps in parallel, in order to assess the impact of specific choices of implementation and modelling. To be considered robust, no results should depend strongly on the specifics of a given componentseparation algorithm. Considering specific details, we generally consider Commander to be the preferred solution on large and intermediate angular scales, due to its somewhat lower largescale effective polarization noise (Table 1), weaker crosscorrelation with the highpass filtered WMAP K−Ka band synchrotron template (Table 2), lower Npoint correlation function fluctuations (Fig. 18 and Table 4), and weaker cosmological parameter dependence on ℓ_{max}, suggesting less internal tension between low, intermediate, and high multipoles (Fig. 15). In addition, the method is able to propagate uncertainties from the input maps to final products by means of Monte Carlo samples drawn from the full posterior. For these reasons, we adopt the Commander solution for the lowℓPlanck 2015 temperature likelihood (Planck Collaboration XI 2016).
However, at high multipoles the Commander solution exhibits a significantly higher effective point source amplitude than the other three maps, due to the exclusion of frequencies below 217 GHz. For temperature, the lowest residual highℓ foregrounds are instead obtained by SMICA, as shown in Fig. 14. As a result, as in 2013, we confirm our preference for the SMICA map for analyses that require fullresolution observations in temperature, such as f_{NL} (first three rows in Table 5 and Planck Collaboration XVII 2016) or lensing reconstruction (Fig. 19 and Planck Collaboration XV 2016). SEVEM is also a very good choice for temperature, since it provides the map with the lowest level of noise at a wide range of scales, as well as a smooth noise power spectrum as measured by the HMHD maps (see Table 1 and Fig. 14). It also performs equally well as SMICA with regard to the estimation of f_{NL}.
In polarization, NILC and SMICA perform equivalently at high multipoles (Fig. 14). The NILC polarization maps yield measurements of f_{NL} that are most consistent with zero (Table 5). The NILC and SMICA analyses also provide an effective mapping of the weights of the Planck frequencies in the needlet/harmonic domains: a given weight tends to be high if a given frequency channel, in a given band in the harmonic domain, is relevant for foreground cleaning; on the other hand, the higher the statistical noise at a given frequency, the lower the associated weight. For NILC and SMICA, the weights for T, E, and Bmodes are shown in Figs. B.2 and D.1, respectively. SEVEM provides the most stability with respect to the effect of bright point sources in polarization (see Sect. 7.1). This is due to the inpainting procedure applied to these sources, which significantly reduces the effect of this contaminant. Therefore, SEVEM could be a suitable choice for those analyses in polarization that cannot easily deal with the presence of point source holes in a mask.
Finally, we note that the SEVEM approach is unique in its ability to provide independent CMB estimates in a number of frequency channels. For analyses that benefit significantly from, or even require, such information, SEVEM is the only meaningful choice. Specific examples include various isotropy estimators (Planck Collaboration XVI 2016), the integrated SachsWolfe stacking analysis (Planck Collaboration XXI 2016), Doppler boosting (Planck Collaboration XXVII 2014), and Rayleigh scattering analyses (Lewis 2013).
10. Conclusions
We have presented four different foregroundreduced CMB maps in both temperature and polarization, derived from the Planck 2015 observations. These maps are based on the full Planck data, including a total of 50 months of LFI observations and 29 months of HFI observations. The temperature component of these maps represents the most accurate description of the CMB intensity sky published to date. In the polarization component, the characteristic Emode signal expected in a standard ΛCDM model is easily discernible over the full sky. Corresponding astrophysical foreground products are described in Planck Collaboration X (2016).
The CMB maps presented here are the direct result of the detailed analyses of systematic errors described in Planck Collaboration VI (2016) and Planck Collaboration VIII (2016), which led to an effective reduction of systematic errors by almost two orders of magnitude in power on large angular scales in polarization compared to 2013. However, despite these improvements, the polarization systematic errors in the Planck 2015 data set are not yet negligible in several frequencies on the very largest scales. Multipoles below ℓ = 20 are therefore suppressed by lowpass filtering in the current componentseparated CMB maps.
Additionally, as already noted in Planck Collaboration XII (2016), we observe a mismatch in the effective noise amplitude of 10–20% when comparing the latest generation of Planck simulations (FFP8) with the data. Considering that the current temperature sky maps are signaldominated up to ℓ ≈ 2000, this noise mismatch is of little practical importance for cosmological analyses based on temperature observations, except on the very smallest scales. As in 2013, the temperature maps presented here are therefore used for a wide range of important applications, including largescale temperature likelihood estimation (Planck Collaboration XI 2016), gravitational lensing (Planck Collaboration XV 2016), studies of isotropy and statistics (Planck Collaboration XVI 2016), primordial nonGaussianity (Planck Collaboration XVII 2016), and nontrivial cosmological topologies (Planck Collaboration XVIII 2016). For highℓ power spectrum and likelihood estimation, we recommend the crossspectrum based methods described in Planck Collaboration XI (2016), primarily due to difficulties in establishing sufficiently accurate models of unresolved extragalactic highℓ foregrounds for the maps presented here. Cosmological parameters derived by temperature power spectra from these maps have been compared with the results of the Planck 2015 likelihood analysis (Planck Collaboration XI 2016) and found to agree at the 1σ level.
For polarization, the noise mismatch is not negligible, and we therefore do not yet recommend using the provided maps for cosmological studies that require a highly accurate noise model. The polarization maps presented here may still be very useful for many important cosmological applications, including crosscorrelation and crossspectrum based analyses, and we therefore release the maps despite the current noise mismatch. Analysis of the higherorder statistics of these maps has been performed within the current framework of precision assessment and presented in this paper. The bispectrum analysis, including Emode polarization, applied to the results of the four component separation methods, gives evidence of a vanishing nonGaussian signal for three geometrical configurations, namely local, orthogonal, and equilateral. The Bmodes derived from the current polarization maps have been crosscorrelated with the predicted lensing Bmodes from the measured E signal and the lensing potential measured independently from temperature in Planck Collaboration XV (2016); the result is found to be in excellent agreement among the four component separation methods, as well as with the prediction of the Planck bestfit cosmology.
On the basis of these encouraging results, intense work is continuing to reduce the largescale polarization systematics to negligible levels, as well as to resolve the noise simulation mismatch, and good progress is being made. Updated products will be published as soon as this work has reached a successful completion.
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).
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, 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/planckcollaboration. Some of the results in this paper have been derived using the HEALPix package.
References
 Basak, S., & Delabrouille, J. 2012, MNRAS, 419, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Basak, S., & Delabrouille, J. 2013, MNRAS, 435, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Cardoso, J., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 735 [Google Scholar]
 Cruz, M., Vielva, P., MartínezGonzález, E., & Barreiro, R. B. 2011, MNRAS, 412, 2383 [NASA ADS] [CrossRef] [Google Scholar]
 Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2005, ApJ, 622, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Dickinson, C., Lawrence, C. R., et al. 2006, ApJ, 641, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., O’Dwyer, I. J., Jewell, J. B., et al. 2004, ApJS, 155, 227 [NASA ADS] [CrossRef] [Google Scholar]
 FernándezCobos, R., Vielva, P., Barreiro, R. B., & MartínezGonzález, E. 2012, MNRAS, 420, 2162 [NASA ADS] [CrossRef] [Google Scholar]
 Finkbeiner, D. P. 2003, ApJS, 146, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Gjerløw, E., Eriksen, H. K., Banday, A. J., Górski, K. M., & Lilje, P. B. 2010, ApJ, 710, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Gold, B., Odegard, N., Weiland, J. L., et al. 2011, ApJS, 192, 15 [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Haslam, C. G. T., Salter, C. J., Stoffel, H., & Wilson, W. E. 1982, A&AS, 47, 1 [NASA ADS] [Google Scholar]
 Komatsu, E., Spergel, D. N., & Wandelt, B. D. 2005, ApJ, 634, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, D., Weiland, J. L., Hinshaw, G., & Bennett, C. L. 2015, ApJ, 801, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Leach, S. M., Cardoso, J.F., Baccigalupi, C., et al. 2008, A&A, 491, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leahy, J. P., Bersanelli, M., D’Arcangelo, O., et al. 2010, A&A, 520, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lewis, A. 2013, J. Cosmol. Astropart. Phys., 8, 53 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., & Bridle, S. 2002, Phys. Rev. D, 66, 103511 [NASA ADS] [CrossRef] [Google Scholar]
 Monteserín, C., Barreiro, R. B., Vielva, P., et al. 2008, MNRAS, 387, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Narcowich, F. J., Petrushev, P., & Ward, J. D. 2006, SIAM J. Math. Anal., 38, 574 [CrossRef] [Google Scholar]
 Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration ES. 2015, The Explanatory Supplement to the Planck 2015 results, http://wiki.cosmos.esa.int/planckpla/index.php/Main_Page (ESA) [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2016, A&A, 594, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2016, A&A, 594, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2016, A&A, 594, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2016, A&A, 594, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2016, A&A, 594, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2016, A&A, 594, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2016, A&A, 594, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2016, A&A, 594, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2016, A&A, 594, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2016, A&A, 594, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2016, A&A, 594, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2016, A&A, 594, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2016, A&A, 594, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2016, A&A, 594, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2016, A&A, 594, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2016, A&A, 594, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2016, A&A, 594, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2016, A&A, 594, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2016, A&A, 594, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2016, A&A, 594, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XIV. 2014, A&A, 564, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rocha, G., Contaldi, C. R., Bond, J. R., & Górski, K. M. 2011, MNRAS, 414, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Rocha, G., Contaldi, C. R., Colombo, L. P. L., et al. 2010, ArXiv eprints [arXiv:1008.4948] [Google Scholar]
Appendix A: Bayesian parametric fitting
Commander (Eriksen et al. 2004, 2008) fits a physical model to a set of observations within a standard Bayesian parametric framework, defined by a set of explicit physical parameters and priors. The code can be operated in two modes, either employing Gibbs sampling to map out the full parameter posterior, or using iterative nonlinear searches to derive the maximumlikelihood solution; the only implementational difference between the two is whether to sample from or maximize the conditional posterior in each Gibbs step. All maps presented in this paper are derived in the maximumlikelihood mode, and uncertainties are evaluated through simulations.
Commander forms the core of the Planck 2015 foregroundtargeted diffuse component separation efforts, and the corresponding results are described in full detail in Planck Collaboration X (2016). In this section we only summarize the most relevant steps for CMBoriented analysis.
Appendix A.1: Intensity
Fig. A.1 Commander processing masks for temperature (top) and polarization (bottom). For temperature, the different shades of grey correspond to different angular resolutions, ranging from 5′ (light grey) through 7.5′ (dark grey) to 40′ FWHM (black). For polarization, the same mask is used for both 10′ and 40′ FWHM resolution. 
In 2013, the Commander CMB temperature solution was derived using only the seven lowest Planck frequency maps between 30 and 353 GHz, adopting a simple fourcomponent signal model, including CMB, CO, modified blackbody thermal dust, and a single powerlaw lowfrequency component (Planck Collaboration XII 2014). The only instrumental parameters included in the analysis were monopole and dipole corrections. In the current release, significantly more data are included in the analysis, and the astrophysical and instrument models have been expanded to account for more effects. Specifically, a total of 32 maps are considered in the analysis, including 21 Planck detector and detectorset maps, 10 WMAP differencing assembly maps, and a 408 MHz lowfrequency survey map (Haslam et al. 1982). This wide frequency range allows us to fit separately for synchrotron, freefree, spinning dust, thermal dust, and CMB, as well as individual CO transitions at 115, 230, and 345 GHz, a common line emission component in the 94 and 100 GHz channels (primarily HCN), and thermal SunyaevZeldovich emission for the Coma and Virgo clusters. On the instrumental side, we now fit for both calibration and bandpass uncertainties, in addition to monopoles and dipoles (Planck Collaboration X 2016).
Fig. A.2 Multipole moment weights used for multiresolution hybridization in the Commander CMB map, as described by Eq. (A.1). 
Fig. A.3 5° × 5° zoomin of the multiresolution contributions to the Commander hybrid CMB map from the 40′ (top left), (top right), and approximately 5′ (bottom left) solutions, centred on the south galactic pole. The hybrid map is shown in the bottom right panel. 
However, the current Commander implementation requires uniform angular resolution across frequencies in order to estimate spectral parameters correctly. This implies that all frequency maps must be smoothed to the resolution of the lowest resolution channel before analysis. In the 2013 Planck release, this problem was partially solved by first determining spectral parameters at low resolution using Commander, as described above, and then solving for the component amplitudes from full resolution data using a socalled Ruler step, leading to the CommanderRuler hybrid. In the current release, we adopt a Commanderonly multistage approach, in which the system is solved at four different angular resolutions, using different subsets of the data, but each with internally coherent angular resolutions. At each higherresolution step, all frequency channels with lower angular resolution are dropped from the analysis, and the model is simplified as required by the new frequency subset. We consider the following four cases in our analysis:

1.
Planck+WMAP+408 MHz at1°FWHM resolution, with a full foreground and instrumental model as described above. A corresponding mask is defined through a χ^{2} evaluation, as described in Planck Collaboration X (2016), and allows 93% of the sky for analysis. This solution represents the most complete astrophysical model of the microwave frequencies on large angular scales available as of today, and forms the basis of the Planck 2015 temperature likelihood at large angular scales (Planck Collaboration XI 2016).

2.
Planck 30–857 GHz at 40′FWHM resolution, with a single powerlaw foreground model at low frequencies to account simultaneously for synchrotron, freefree and spinning dust emission. All global instrumental parameters are fixed to the values obtained in the first step. As in the first step, the mask is defined in terms of the perpixel χ^{2}, but with a higher (i.e., less restrictive) threshold, resulting in a total sky fraction of 98.4%; see Fig. A.1. The lowfrequency foreground sector of this model is identical to that adopted for the Planck 2013 component separation analysis.

3.
Planck 143–857 GHz at FWHM resolution, including only CMB, CO, and thermal dust emission. The dust emissivity is fitted in each pixel, while the dust temperature is fixed to the lowresolution solution. The mask is defined as the product of the perpixel χ^{2} at the native resolution, the 40′ FWHM mask, and an external point source mask, and it admits a sky fraction of 95%.

4.
Planck 217–857 GHz at ~ 5′FWHM resolution, fitting only the CMB and thermal dust emission amplitude per pixel. All other parameters are fixed to their lowerresolution values. The highestresolution mask is defined in the same way as the FWHM one, but with an additional cut on the CO amplitude of 0.5 K km s^{1} (see Planck Collaboration X 2016), since CO is no longer part of the foreground model. This mask admits a sky fraction of 82.4%.
A comparison of the various Planckonly masks are shown in Fig. A.1, while the likelihood mask is shown in Planck Collaboration X (2016).
Fig. A.4 Difference between the Commander CMB solutions derived using (1) only Planck data and (2) Planck+WMAP+408 MHz, both smoothed to a common resolution of 1° FWHM. The grey region indicates the common mask. 
Given these singleresolution maps, we hybridize the three Planckonly solutions (ranging between 40 and 5′ FWHM) into a single map as follows. First, we replace the masked regions at each resolution with a constrained Gaussian realization (see Eriksen et al. 2004 for details), drawn from P(a_{CMB}  C_{ℓ},d), to suppress ringing and Galactic leakage effects. Then, we combine the three Planckonly maps in harmonic space using the following cosine apodization scheme, (A.1)where is the beam transfer function and is the pixel window function for a particular resolution, and the cosine apodization weights are given by as illustrated in Fig. A.2. Figure A.3 shows the contributions from the three different resolutions and their sum.
The result is a single fullresolution CMB map with variable effective sky fraction as a function of multipole. For multipoles below ℓ = 200, a total of 98.4% of the sky are formally derived from real measurements, while the remaining 1.6% are filled with a Gaussian constrained realization. For multipoles between ℓ = 300 and 1000, 95% of the sky are derived from real data, while at multipoles above ℓ = 1200, only 82% of the sky are derived from real data. However, the fact that the map is datadriven in a particular pixel for a given angular scale does not necessarily guarantee proper goodnessoffit, and we therefore recommend application of the most conservative 82% mask for highprecision science. For analyses requiring maximum sky coverage, we instead recommend the likelihood map, which supports 93% of the sky, albeit at lower angular resolution. Figure A.4 shows the difference between the Planckonly and the Planck+WMAP+408 MHz Commander CMB solutions at 1° FWHM resolution. The rms difference between these two maps is 2.7 μK outside the common UT78 mask.
Appendix A.2: Polarization
The Commander CMB polarization map is derived in an analogous manner to the temperature solution, but relying on Planck observations alone. At low resolution, we derive a 40′ map from all frequencies between 30 and 353 GHz, including CMB, synchrotron and thermal dust in the signal model. At high resolution, we derive a 10′ map from frequencies between 100 and 353 GHz, including only CMB and thermal dust. The two maps are hybridized into one 10′ map using the same formalism as in Eq. (A.1), but with only one hybridization operator covering ℓ = 200−300.
As discussed in Planck Collaboration X (2016), several spectral index models have been explored for polarized synchrotron and thermal dust emission, including: (1) spatially constant; (2) smoothed over large angular scales; (3) based on the temperature model. The main conclusion, however, is that the current polarization data have too low signaltonoise ratio for both synchrotron and thermal dust indices to discriminate between the three models at a statistically significant level. Allowing such additional degrees of freedom only increases the degeneracies between the various components without improving the overall fit, and also leaves the solution sensitive to largescale residual systematics. For now, we therefore adopt the temperaturederived spectral parameters for both synchrotron and thermal dust for the primary Commander CMB polarization map.
For polarization masking purposes, we construct a mask from the product of the thresholded lowresolution χ^{2} map and the CO amplitude map. The same mask is applied at both the 40 and 10′ FWHM resolutions.
Appendix B: Internal linear combination in needlet space
The goal of NILC is to estimate the CMB from multifrequency observations while minimizing the contamination from Galactic and extragalactic foregrounds, as well as instrumental noise. The method makes a linear combination of the data from the input maps with minimum variance using a basis of spherical wavelets called needlets (Narcowich et al. 2006). Due to their unique properties, needlets allow localized filtering in both pixel space and harmonic space. Localization in pixel space allows the weights of the linear combination to adapt to local conditions of foreground contamination and noise, whereas localization in harmonic space allows the method to favour foreground rejection on large scales and noise rejection on small scales. Needlets permit the weights to vary smoothly on large scales and rapidly on small scales, which is not possible by cutting the sky in zones prior to processing.
The NILC pipeline (Basak & Delabrouille 2012, 2013) is applicable to scalar fields on the sphere, hence we work separately on maps of temperature and the E and Bmodes of polarization. The decomposition of input polarization maps into E and B is done on the full sky. At the end, the CMB Q and U maps are reconstructed from the E and B maps.
Prior to applying NILC, all of the input maps are convolved or deconvolved in harmonic space to a common resolution corresponding to a Gaussian beam of 5′ FWHM. Each map is then decomposed into a set of needlet coefficients. For each scale j, needlet coefficients of a given map are stored in the form of a single HEALPix map. The filters used to compute filtered maps are shaped as follows For each scale j, the filter has compact support between the multipoles and with a peak at (see Table B.1 and Fig. B.1). The needlet coefficients are computed from these filtered maps on HEALPix pixels with N_{side} equal to the smallest power of 2 larger than .
List of needlet bands used in the NILC analysis.
Fig. B.1 Needlet bands used in the analysis. The thick black line shows the normalization of the needlet bands, that is, the total filter applied to the original map after needlet decomposition and synthesis of the output map from needlet coefficients. 
Fig. B.2 Fullsky average of needlet weights for different frequency channels and needlet bands. From top to bottom, the panels show results for temperature, Emodes, and Bmodes. 
In order to show the contribution of the various frequency channels to the final CMB map at different needlet bands, we compute the full sky average of needlet weights for each frequency channel and needlet band. Figure B.2 shows that most of the contribution to the reconstructed CMB maps comes from the 143 GHz and 217 GHz channels. In the lowℓ needlet bands, the contribution from 143 GHz is large compared to that from 217 GHz. However, due to better angular resolution, the 217 GHz channel contributes more than the 143 GHz channel in the highest ℓ needlet bands. In the intermediate needlet bands, the contributions from these two channels are comparable. It is interesting to note the contribution from the LFI in the lowest ℓ bands; in intensity, the 70 GHz channel serves mostly for foreground removal, while for polarization it contributes to the CMB solution with a positive weight, and very similarly between E and B. We stress again that the lowest ℓ modes have been filtered out in the results presented here, and therefore the low ℓ results will require further investigation.
Appendix B.1: Masking
The confidence masks for NILC for intensity and polarization have been generated following a procedure similar to that used by SMICA, but adopting a different parameterization.
For intensity, the NILC CMB map is filtered through a spectral window (B.1)The result is then squared and smoothed with a Gaussian circular beam with FWHM 120′. The variance map obtained in this way is then corrected for the noise contribution by subtracting the variance map for the noise obtained in the same way from the NILC HRHD map. The confidence map is obtained by thresholding the noisecorrected variance map at 73.5 μK^{2}. For polarization, the polarized intensity is obtained from the NILC outputs, and is filtered through a circularly symmetric Gaussian window function of FWHM 30′. The result is then squared and smoothed with a Gaussian of FWHM 210′. The resulting variance map is corrected for the noise contribution by following the same procedure used for intensity. The confidence map is obtained by thresholding at 6.75 μK^{2}. The resulting masks are shown in Fig. B.3.
Fig. B.3 NILC masks for temperature (top) and polarization (bottom). 
Appendix C: Template fitting
The SEVEM method (Leach et al. 2008; FernándezCobos et al. 2012) aims to produce clean CMB maps for several frequency channels by using internal template fitting. The templates are constructed from the Planck data, typically as the subtraction of two close Planck frequency channels to remove the CMB signal. The beam sizes of the maps are equalized before subtraction. The cleaning is achieved simply by subtracting a linear combination of the templates t_{j}(x) from the data, with coefficients α_{j} obtained by minimizing the variance outside a given mask, (C.1)Here n_{t} is the number of templates used, and T_{c}(x,ν) and d(x,ν) correspond to the cleaned and raw maps at frequency ν, respectively. The same expression applies for T, Q, or U.
The cleaned frequency maps are then combined in harmonic space, taking into account the noise level, resolution, and, for temperature, a rough estimate of the foreground residuals of each cleaned channel, to produce a final CMB map at the required resolution.
Appendix C.1: Implementation for temperature
For temperature, we followed a similar procedure to that used for the Planck 2013 release, i.e., the 100, 143, and 217 GHz maps are cleaned using four templates constructed from the six remaining frequency channels. A few differences have been implemented in the current pipeline with respect to the previous work: the use of a single coefficient over the whole sky for each template (instead of defining two regions), the use of inpainting to reduce contamination from sources, and the use of the 857 GHz channel as a template (instead of 857−545). We note that the other three templates (30−44, 44−70, and 353−217) are the same as for the previous release.
The six frequency channels used to construct templates (30 to 70 GHz and 353 to 857 GHz) are inpainted at the position of sources detected by the Mexican hat wavelet (MHW) algorithm (Planck Collaboration XXVI 2016). The size of the holes to be inpainted is determined by taking into account the beam size of the channel as well as the flux density of each source. We do a simple diffusive inpainting, which fills one pixel with the mean value of the neighbouring pixels in an iterative way. To avoid inconsistencies when subtracting two channels, each map is inpainted on the sources detected in both that map and on the other map used to construct the template. For example, to construct the (30−44) template, both maps are inpainted in the positions of the sources detected at 30 and 44 GHz. This reduces significantly the contamination from compact sources in the templates.
Once they have been inpainted, the maps are smoothed to a common resolution and then subtracted. To construct the first three templates, the first channel in the subtraction is smoothed with the beam of the second map and vice versa. For the 857 GHz template, we simply filter the map with the 545 GHz beam (this is done for comparison with the 857−545 template from the 2013 pipeline, where the 857 GHz map was smoothed by the 545 GHz beam). We note that the coefficients used to multiply this template are typically of order 10^{4}, so the level of CMB signal introduced by this template in the final cleaned map is negligible. We take advantage of this to drop subtraction of the 545 GHz map, as was done in 2013, nominally to remove the CMB signal. This simplifies the method and also reduces the noise in this template.
The 100, 143, and 217 GHz maps are then cleaned by subtracting a linear combination of the four templates. The coefficients of the linear combination (Table C.1) are obtained by minimizing the variance outside an analysis mask. The main difference with respect to the 2013 release is that we have used the same coefficients for the whole sky (instead of dividing it into two regions), since this simplifies the procedure without affecting the quality of the reconstruction (other than on those pixels very close to the Galactic centre that need to be masked in any case). This analysis mask covers the 1% of the sky with the brightest emission, as well as sources detected at all frequency channels. Once the maps are cleaned, each is inpainted on the source positions detected at that (raw) channel. Then, the MHW algorithm is run again, now on the cleaned maps. A relatively small number of new sources are found and are also inpainted at each channel. The resolution of the cleaned map is the same as that of the raw map. We note that no assumptions about the noise or foregrounds are made in order to construct the singlefrequency cleaned maps.
Linear coefficients, α_{j}, of the templates used to clean individual frequency maps with SEVEM for temperature.
Linear coefficients, α_{j}, of the templates used to clean individual frequency maps with SEVEM for Q and U.
Finally, the SEVEM CMB map is constructed by combining the cleaned and inpainted 143 and 217 GHz maps. In the combination, the maps are weighted taking into account the noise, resolution, and a rough estimate (obtained from realistic simulations) of the foreground residuals in each map. The weights are shown in Fig. C.1. The resolution of this map corresponds to a Gaussian beam of FWHM 5′ and HEALPix resolution N_{side} = 2048, with a maximum multipole ℓ_{max} = 4000.
Fig. C.1 Weights used to combine the 143 and 217 GHz cleaned frequency channel maps into the final SEVEM temperature map. The weights do not sum to unity because they include the effect of deconvolving the beams of the frequency maps and convolving with the 5′ beam of the final map. 
The same procedure, including the full inpainting process for point sources, is applied when running the pipeline on FFP8 simulations.
Appendix C.2: Implementation for polarization
To clean the polarization maps, a procedure similar to the one used for the temperature maps is applied to the Q and U maps independently. However, given that narrower frequency coverage is available for polarization, the templates and maps to be cleaned are different. In particular, we clean the 70, 100, and 143 GHz maps using three templates. Since the signaltonoise ratio is lower, at 100 and 143 GHz, the clean maps are produced at N_{side} = 1024, with resolution corresponding to a Gaussian beam of FWHM 10′ and a maximum multipole ℓ_{max} = 3071. At 70 GHz, the map is produced at its native resolution.
The first step of the pipeline is to inpaint the positions of the sources detected using the MHW algorithm in those channels that will be used to construct templates. Similarly to the temperature case, for a given template (constructed as the difference of two frequency channels), the inpainting is performed for all of the sources detected in any of the channels. We note that the inpainting is performed in the frequency maps at their native resolution.
These inpainted maps are then used to construct a total of four templates. To trace the synchrotron emission, we construct a 30−44 template. For the dust emission, the following templates are used: 353−217 (smoothed to 10′ resolution); 217−143 (used to clean 70 and 100 GHz); and 217−100 (to clean 143 GHz). These last two templates are constructed at 1° resolution, since an additional smoothing becomes necessary in order to increase the signaltonoise ratio of the template. Otherwise, the estimated coefficients are driven by the noise and the cleaning is less efficient. Since fewer frequency channels are available in polarization, it becomes necessary to use the maps that are to be cleaned as part of one of the templates. Therefore the 100 GHz map is used to clean the 143 GHz frequency channel and vice versa, making the clean maps less independent than in the temperature case.
These templates are then used to clean the noninpainted 70 (at its native resolution), 100 (at 10′ resolution), and 143 GHz maps (also at 10′). The corresponding linear coefficients (listed in Table C.2) are estimated independently for Q and U by minimizing the variance of the clean maps outside a mask that covers compact sources and the 3% of sky with the brightest Galactic emission. Once the maps have been cleaned, inpainting of the sources detected at each map is carried out. The size of the holes to be inpainted takes into account the additional smoothing of the 100 and 143 GHz maps. Similarly to temperature, the same inpainting processing is applied to the point source positions when running the pipeline on the FFP8 simulations. The 100 and 143 GHz clean maps are then combined in harmonic space, using a fullsky E and B decomposition, to produce the final CMB maps for the Q and U components with a Gaussian beam of FWHM 10′ and HEALPix resolution N_{side} = 1024. Each map is weighted taking into account its noise level at each multipole.
Before applying the postprocessing highpass filter to the cleaned Q and U polarization maps, we inpaint the region with the brightest Galactic residuals (5% of the sky) with the same simple algorithm used for point source holes. This is done to avoid introducing ringing around the Galactic centre when the maps are filtered.
In addition, E and Bmode maps are constructed from the clean Q and U maps. In this case, prior to performing the decomposition, the region of the Q and U maps defined by the SEVEM confidence mask is filled with a Gaussianconstrained realization (Eriksen et al. 2004).
Appendix C.3: Masks
Fig. C.2 SEVEM masks in temperature (top) and polarization (bottom). 
In temperature, the SEVEM confidence mask is produced by thresholding differences obtained between different CMB reconstructions. In particular, we construct the difference map between the clean 217 and 143 GHz maps at a resolution of FWHM 30′ and N_{side} = 256. The brightest pixels of this map (and their direct neighbours) are successively removed from the CMB combined map (at the same resolution) and the dispersion of the CMB combined map calculated. If a sufficient number of pixels is removed, the dispersion of the CMB map goes down and reaches a plateau, indicating that convergence has been achieved. The removed pixels constitute the mask. This mask is then smoothed with a Gaussian beam of 1° to avoid sharp edges, and upgraded to full resolution. The same procedure is repeated for the other two difference maps: the cleaned 143−100 map; and the difference of two clean CMB combined maps, whose linear coefficients have been obtained by minimizing the variance outside two different masks. Finally, the three masks are multiplied in order to produce our final confidence mask, which leaves a suitable sky fraction of approximately 85%. Residual monopoles and dipoles outside this mask are subtracted for the single frequency and combined cleaned maps.
For polarization, the clean combined map is downgraded to a resolution of FWHM 90′ and N_{side} = 128. The dispersion at each pixel is estimated from a circle centred in the pixel being considered. Those pixels with a dispersion above a given threshold are included in the mask, which is then smoothed with a Gaussian beam of FWHM 90′ to avoid sharp edges and upgraded to the required N_{side}. Finally, this mask is multiplied by a mask customized to cover the CO emission, in order to discard those pixels contaminated by this foreground component due to the bandpass leakage. An additional 1% of the sky is added to the final mask to remove those pixels most affected by the highpass filtering that is subsequently applied to the cleaned maps. The final mask allows for a useful fraction of the sky of approximately 80%. The 2015 SEVEM masks are shown in Fig. C.2.
Appendix D: Spectral matching
Spectral Matching Independent Component Analysis (SMICA) is a semiblind component separation algorithm that operates in harmonic space. CMB maps are synthesized from spherical harmonic coefficients, s_{ℓm}, obtained as weighted linear combinations of the coefficients of N_{chan} input maps, (D.1)where x_{ℓm} is the N_{chan} × 1 vector of the spherical harmonic coefficients of the input maps and w_{ℓ} is the N_{chan} × 1 vector of weights. Here and in the following, for the sake of exposition, a realvalued spherical harmonic basis is used (it is always possible to do so) so that all the quantities are real. In practice, we use the standard complex basis, but this is immaterial. The spectral weights used to produce the SMICA maps are designed to minimize the total foreground and noise contamination at each multipole, under the constraint that the resulting map has a well defined effective beam window function, that of a Gaussian beam with 5′ FWHM.
In theory, such weights are given by (D.2)where the N_{chan} × 1 vector a is the frequency spectrum of the CMB and the N_{chan} × N_{chan} spectral covariance matrix R_{ℓ} contains in its (i,j) entry the crosspower spectrum at multipole ℓ between the input frequency maps i and j, at 5′ resolution. In practice, the spectral covariance matrices must be estimated from the data. The natural sample estimate (D.3)can be used, possibly after some binning, to replace R_{ℓ} in the weight formula, Eq. (D.2). This works well at high ℓ because a large number of modes are averaged in Eq. (D.3), so this estimate has low variance. At larger scales, however, it is necessary to constrain the spectral covariance matrices in order to obtain reliable estimates and mitigate the effects of chance correlations. For that purpose, SMICA uses a semiparametric model of these matrices.
The CMB, the foregrounds, and the noise are independent processes, so the spectral covariance matrices, after beam correction, can be decomposed into (D.4)We assume that the noise is uncorrelated between frequency channels, therefore the noise contribution is diagonal. The signal parts are modelled as (D.5)where C_{ℓ} is the CMB power spectrum, matrix F is N_{chan} × d, and P_{ℓ} is a d × d positive definite matrix. A SMICA fit consists of fitting the model of Eqs. (D.4) and (D.5) to empirical covariance matrices by minimizing the spectral matching criterion The fitted values of R_{ℓ} can be seen as regularized versions of their empirical counterparts , to be used in Eq. (D.2). If no constraints are imposed on matrices F and P_{ℓ}, except that the latter is positive definite, the foreground model is equivalent to d sky templates with arbitrary frequency spectra (represented by the columns of F) and arbitrary power spectra and correlations between them (represented by P_{ℓ}). Ultimately, the foreground contribution is controlled by a single parameter, the rank d of the foreground model.
More details regarding the principles of SMICA can be found in Cardoso et al. (2008). We now describe its extension to polarization data.
There are several options to extend SMICA to polarization. We choose to obtain the CMB Q and U maps through a joint processing of the E and Bmodes. The CMB U and Q maps are synthesized from their E and Bmodes, and , obtained as linear combinations (D.6)where the N_{chan} × 1 vectors and are the spherical harmonic coefficients of the input maps and W_{ℓ} is a 2N_{chan} × 2 matrix of weights. The optimal weights are obtained as a simple generalization of the singlefield case, (D.7)where R_{ℓ} now is a 2N_{chan} × 2N_{chan} covariance matrix. The weights defined by Eq. (D.7) can be safely obtained at high ℓ by using the sample covariance matrices, possibly after binning them. At low multipoles, some regularization via modelling is in order, as for the temperature analysis. We use a natural extension of temperature Eq. (D.5), in the form where F_{E} and F_{B} are N_{chan} × d matrices and P_{ℓ} is a 2d × 2d matrix. For this release, no constraints are imposed on those matrices, except that each P_{ℓ} must be positive definite. Hence, as in temperature, SMICA fits a foreground model representing d polarized templates with arbitrary frequency spectra, power spectra, and correlations.
Fig. D.1 SMICA weights for temperature (top) and polarization (bottom). For readability, the values are shown for input maps in units of antenna temperature. The plot goes up to ℓ ≃ 3600, but the output maps are synthesized using all multipoles up to ℓ = 4000. For polarization, the thick solid lines show the contribution of input Emodes to the CMB Emodes and the thick dashed lines show the same for the Bmodes. The thin lines, all close to zero, show “crosscontributions” of input Emodes to the CMB Bmodes and vice versa. 
Appendix D.1: Implementation for temperature
The production of the CMB temperature map is mostly unchanged from 2013 (see Planck Collaboration XII 2014 for details). Here, we recall the 3step fitting strategy adopted in 2013 and mention the changes made for this release.

First fit: recalibration. A preliminary and independent SMICA fit is done with power spectra estimated over a clean fraction of the sky, including a as a free parameter. This can be understood as a recalibration procedure; the estimated value of a is kept fixed in the following steps. The number of foreground templates, d, is now 5, whereas it was 4 in 2013. The 30 GHz channel is not recalibrated, unlike in 2013.

Second fit: foreground emission. The foreground emission, matrix F, is estimated in a second SMICA fit with a kept fixed at the value found in the first step. Spectra are estimated over a large fraction of the sky (f_{sky} = 97%) and the fit is made over the range 4 ≤ ℓ ≤ 150. For the data, this step is the same as in 2013. For FFP8, foreground emission is captured using seven templates, whereas six templates were used for the FFP6 simulations in 2013.

Third fit: power spectra. The frequency spectra captured by vector a and matrix F are kept fixed at the values found in steps 1 and 2; we fit only the signal power spectra C_{ℓ} and P_{ℓ} and the frequency channel noise power spectra. Spectral covariance matrices are computed over 97% of the sky and the fit includes all multipoles up to ℓ = 1500. This step is the same as in 2013.
Fig. D.2 SMICA masks in temperature (top) and polarization (bottom). 
Appendix D.2: Implementation for polarization
All Planck polarized channels are used to produce the polarized CMB maps. The process is easier in polarization than it is in temperature because less precision is required due to the lower signaltonoise ratio and also because the foregrounds appear to have a simpler structure. In particular, we do not preprocess the frequency maps by subtracting or masking point sources, as is done to the temperature maps.
The SMICA fit in polarization is conducted with the same parameters as in temperature, but with two differences. First, the recalibration step is omitted, and instead we use the CMB frequency spectrum (vector a) determined from temperature maps. Second, the foreground model comprises six polarized templates (as for temperature) but the matrices F_{E} and F_{B} are fitted in the second step over the multipole range 4 ≤ ℓ ≤ 50. The weights determined for the polarization data are shown at the bottom of Fig. D.1.
In order to mitigate spectral leakage from E to B and from the Galactic plane onto other regions of the sky, we do not compute spherical transforms directly from masked Q and U maps. Instead, for each pair of input Q and U maps, we first produce fullsky E and B maps to which an apodized Galactic mask is applied. It is from those masked E and B maps that spherical harmonic coefficients are computed for the estimation of the spectra and crossspectra going into and for the synthesis of the final CMB map.
Appendix D.3: Masks
Confidence masks, shown in Fig. D.2, are built using the procedure described in Planck Collaboration XII (2014), with the following changes. For temperature, we apply a bandpass filter to the CMB map that is then squared and smoothed at 3.̊5 (it was 2° in 2013). The confidence mask is obtained by thresholding the resulting map of local power. The threshold is determined by visual inspection to be 50 μK^{2} (it was 70 μK^{2} in 2013). The resulting mask is then enlarged by multiplication by a Galactic mask covering 10% of the sky. For polarization, the confidence mask is obtained by a procedure similar to temperature, but the CMB Q and U maps are lowpass filtered (rather than bandpassed) using a Gaussian beam with 30′ FWHM. They are then squared and smoothed to 3.̊5 resolution. Any area where the resulting map is above 5 μK^{2} is excluded from the confidence mask. In addition to that, we exclude pixels that are less than 7° away from the Galactic equator.
Appendix E: FFP8 simulations
In this Appendix we provide a compendium of analyses evaluated from the FFP8 data set, corresponding directly to those performed on the Planck 2015 data in the main text.
Appendix E.1: CMB map differences
Figures 4, 9, and 10 show pairwise difference maps between any of the Commander, NILC, SEVEM, and SMICA maps derived from the Planck 2015 data for each of the three Stokes parameters, I, Q, and U. In Figs. E.1–E.3 we show the same quantities evaluated from the FFP8 simulation set.
Fig. E.1 Pairwise differences between CMB temperature maps obtained from FFP8 simulations. Prior to differencing, the maps have been smoothed to 80′ FWHM and downgraded to N_{side} = 128. 
Overall, the relative differences are of similar magnitude in the simulated data set as in the data, although they have somewhat different dominant morphology. The highlatitude differences are in general somewhat weaker in FFP8 than in the data, while the lowlatitude differences are somewhat stronger. This is primarily due to the rather complicated foreground model adopted for thermal dust emission in the FFP8 simulations (see Planck Collaboration XII 2016 for details). In short, thermal dust emission is modelled in the FFP8 simulations by a sum of two modified blackbody components, one of which has a temperaturedependent spectral index, β(T_{d}). This was motivated by the results presented by Planck Collaboration Int. XIV (2014); however, as shown by the updated analysis in Planck Collaboration X (2016), there is no evidence for steepening from the Planck data alone. Only when the IRAS 100μm data are included in the fit is any such effect seen. At such high frequencies the thermal dust physics is far more complicated, and the overall calibration problem more difficult. For the present discussion, it is sufficient to note that the thermal dust model adopted for the FFP8 simulations is significantly more complicated in terms of frequency dependence than the observed sky, and includes significant spectral curvature below 353 GHz that is not seen in the actual data. None of the component separation codes account for this additional curvature, and this results in the strong features near the Galactic plane seen in some of the pairwise difference maps shown in Fig. 4.
Fig. E.2 Pairwise differences between CMB Q maps obtained from FFP8 simulations. Smoothing and degrading is as in Fig. E.1. 
Fig. E.3 Pairwise differences between CMB U maps obtained from FFP8 simulations. Smoothing and degrading is as in Figs. E.1 and E.2. 
Appendix E.2: 1point statistics of the total intensity FFP8 maps
Next, we consider the 1point statistics of the FFP8 fiducial map in temperature, using the FFP8UT74 mask (Table E.1). These statistics were defined in Sect. 7.1, and corresponding polarization results were shown in Table 3. We see that all the methods are in good agreement with the Monte Carlo simulations at all resolutions considered, and the differences among codes are small.
Lower tail probability in percent for the variance, skewness, and kurtosis for the FFP8 total intensity analysis at three different resolutions.
Appendix E.3: The realspace Npoint correlation functions for FFP8 maps
Finally, we present the realspace Npoint correlation functions derived from the FFP8 simulations, complementing the analysis in Sect. 7.2. As we did for the data, we analyse the FFP8 CMB temperature and highpass filtered polarization estimates downgraded to N_{side} = 64 and smoothed with a Gaussian kernel of 160′ FWHM. For the polarization analysis we employ a low resolution version of the common polarization mask, and for the temperature analysis we employ a low resolution version of the common temperature mask.
The Npoint functions for the FFP8 simulations are shown in Figs. E.4, plotted in terms of differences between the Npoint functions for the highpass filtered fiducial FFP8 map and the mean value estimated from 1000 Monte Carlo simulations. The probabilities of obtaining larger values for the χ^{2} statistic, compared to the Planck fiducial ΛCDM model, are tabulated in Table E.2.
Overall, there is good agreement between the recovered FFP8 CMB estimates and the corresponding Monte Carlo ensembles. The main outlier is NILC, for which the fluctuations are larger than expected for all 3point polarization functions.
Appendix E.4: CMB Power spectra and residuals from individual components
Here we discuss the recovery of the power spectrum in the FFP8 simulations, and assess the impact of foregrounds in the CMB solutions presented in this paper. CMB angular power spectra from the foreground cleaned CMB maps from FFP8 are shown in Fig. E.5, compared to the input for TT and EE. The performance on simulations is comparable to that on the data, as may be seen by comparison with Fig. 14.
Figures E.6 and E.7 show the residual effect of combined individual components on the foregroundcleaned CMB maps. These plots give an estimate of how the residuals of the foreground components in the CMB solutions presented in this paper compare in power and are distributed in ℓ. The FFP8 simulated skies of the individual components are propagated through the various componentseparation pipelines, and their power spectra are calculated using the FFP8UT74 sky mask. The individual components are those simulated in FFP8 (Sect. 2.2 and Planck Collaboration XII 2016) and include in TT dust, CIB fluctuations, unresolved sources (PS), and the sum of the remaining components (“other”). In EE, they include synchrotron, dust, noise, and the sum of the remaining components. With the exception of the lowest multipoles, the dominant contaminant to the TT signal is noise. For EE, the residuals from dust dominate over the other foregrounds, but remain subdominant with respect to instrumental noise at all multipoles.
Fig. E.4 Difference between the Npoint functions for the highpass filtered N_{side} = 64 FFP8 CMB estimates and the corresponding means estimated from 1000 MC simulations. The Stokes parameters Q_{r} and U_{r} were locally rotated so that the correlation functions are independent of coordinate frame. The first row shows results for the temperature map, from left to right, results for the 2point, pseudocollapsed 3point, equilateral 3point and connected rhombic 4point functions. We note that the lines lie on top of each other. The second row shows results for the 2point function, from left to right, TQ_{r}, TU_{r}, Q_{r}Q_{r}, and U_{r}U_{r}. The third row shows results for the pseudocollapsed 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}, and the fourth row shows results for the equilateral 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}. The black solid, red tripledotdashed, orange dashed, green dotdashed, and blue long dashed lines correspond to the true, Commander, NILC, SEVEM, and SMICA maps, respectively. The true CMB map was analysed with added noise corresponding to the SMICA component separation method. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using SMICA simulations. See Sect. 7.2 for the definition of the separation angle θ. 
Fig. E.5 Power spectra of the foregroundcleaned CMB maps from FFP8 simulations. Top: TT power spectra evaluated using the FFP8UT74 mask. Bottom: EE power spectra evaluated using the FFP8UP76 mask. Thick lines show the spectra of signal plus noise estimated from the halfmission halfsum maps; thin lines show the noise levels from halfmission halfdifference maps. The black line shows the input spectrum. 
Fig. E.6 TT angular power spectra of residuals from the indicated FFP8 components in the Planck 2015 CMB maps, compared with the predicted signal from the best fit cosmology. “Other” is the sum of CO, freefree, thermal and kinetic SZ, spinning dust, and synchrotron emission. The horizontal axis is linear in ℓ^{ 0.5}. 
Fig. E.7 EE angular power spectra of residuals from the indicated FFP8 components in the Planck 2015 CMB maps, compared with the predicted signal from the best fit cosmology. “Other” is the sum of CO, freefree, thermal and kinetic SZ, spinning dust, CIB fluctuations, and unresolved radio and infrared sources. The horizontal axis is linear in ℓ^{ 0.5}. 
All Tables
Masks and statistics of componentseparated CMB maps from data and FFP8 simulations.
Correlation coefficients of CMB maps with foreground templates for temperature and polarization.
Percentage of simulations showing a lower variance, skewness, or kurtosis than the P maps, the “lowertail probability”.
Linear coefficients, α_{j}, of the templates used to clean individual frequency maps with SEVEM for temperature.
Linear coefficients, α_{j}, of the templates used to clean individual frequency maps with SEVEM for Q and U.
Lower tail probability in percent for the variance, skewness, and kurtosis for the FFP8 total intensity analysis at three different resolutions.
All Figures
Fig. 1 Preferred masks for analysing componentseparated CMB maps in temperature (left) and polarization (right). 

In the text 
Fig. 2 Componentseparated CMB temperature maps at full resolution, FWHM 5′, N_{side} = 2048. 

In the text 
Fig. 3 Differences between the componentseparated CMB temperature maps from the 2013 and the 2015 releases. The maps have been smoothed to FWHM 80′ and downgraded to N_{side} = 128. 

In the text 
Fig. 4 Pairwise difference maps between CMB temperature maps. As in the previous Fig. 3, the maps have been smoothed to FWHM 80′ and downgraded to N_{side} = 128. 

In the text 
Fig. 5 Difference between output and input CMB temperature maps from FFP8 simulations. Smoothing and downgrading is performed as in Fig. 4. 

In the text 
Fig. 6 Componentseparated CMB Q maps at resolution FWHM 10′, N_{side} = 1024. 

In the text 
Fig. 7 Componentseparated CMB U maps at resolution FWHM 10′, N_{side} = 1024. 

In the text 
Fig. 8 20° × 20° patch of the highpass filtered Commander CMB polarization map, centred on the north ecliptic pole, (l,b) = (96°,30°). Each map is pixelized with a HEALPix resolution of N_{side} = 1024, and has an angular resolution of 10′ FWHM. The top row shows Q and U maps derived from the fullmission data set, the middle row shows the corresponding E and B maps, and the bottom row shows the E and Bmodes of the halfring halfdifference (HRHD) map. Note the characteristic “+” and “×” patterns in the Q and U maps, and the clear asymmetry between E and B in the full data set. Also note that the HRHD (null) E map is consistent with both the full and HRHD B maps. 

In the text 
Fig. 9 Pairwise differences between CMB Q maps, after smoothing to FWHM 80′ and downgrading to N_{side} = 128. 

In the text 
Fig. 10 Pairwise differences between CMB U maps, after smoothing and downgrading are performed as in Fig. 9. 

In the text 
Fig. 11 Difference between output and input CMB Q maps from FFP8 simulations. Smoothing and downgrading are performed as in Figs. 9 and 10. 

In the text 
Fig. 12 Difference between output and input CMB U maps from FFP8 simulations. Smoothing and downgrading are performed as in Fig. 11. 

In the text 
Fig. 13 Power spectra of the the foregroundcleaned CMB maps. Left: TT power spectra evaluated using the UT78 mask. Right: EE power spectra evaluated using the UP78 mask. The thick lines show the spectra of the halfmission halfsum maps containing signal and noise. The thin lines show the spectra of the halfmission halfdifference maps, which give an estimate of the noise and some of the residual systematic effects. The black line shows the Planck 2015 bestfit CMB spectrum for comparison. 

In the text 
Fig. 14 CMB TT and EE power spectra for each of the four foregroundcleaned maps. The first panel shows the TT bandpowers after subtracting the bestfit model of residual extragalactic foregrounds. The black points show the frequencymap based Plik spectrum, and the grey line shows the bestfit ΛCDM model from the Planck 2015 likelihood. The second panel shows the residuals of the TT bandpowers after subtracting the bestfit ΛCDM model. The third and fourth panels show the same information for the EE spectra. The bands are the same for all methods, but the points are offset horizontally for clarity. The parameters of the foreground model are marginalized over when estimating parameters (see Fig. 15). 

In the text 
Fig. 15 Comparison of cosmological parameters estimated from the TT and EE spectra computed from the foregroundcleaned CMB maps. The labels on the horizontal axis give the spectrum and the ℓ_{max} used to obtain each set of results. The black points show the parameters obtained from the frequency mapbased Plik likelihood using the same spectrum and ℓ_{max}. For comparison, we also show the values of the parameters obtained with the full Planck 2015 likelihood (including multipoles up to ℓ_{max} = 2500) as the horizontal line surrounded by a grey band giving the uncertainties. The foreground model used for the CMB maps is the methodtailored fullsky model from the FFP8 simulations in TT and a single D_{ℓ} ∝ ℓ^{2} template in EE. 

In the text 
Fig. 16 Polarized intensity variance (left column), skewness (middle column), and kurtosis (right column) evaluated from the FFP8 Monte Carlo simulations (histogram) and from components of the fiducial FFP8 map at N_{side} = 1024 outside the FFP8UPA76 mask. The variance distributions have been normalized to the mean value of the Monte Carlo distributions for visualization purposes. Coloured vertical lines correspond to different combinations of components: the sum of CMB and noise is shown in brown; the sum of CMB, noise, and thermal dust is shown in light teal; the sum of CMB, noise, and radio point sources is shown in violet; the sum of CMB, noise, and all foregrounds is shown in pink. 

In the text 
Fig. 17 Polarized intensity variance evaluated from the FFP8 Monte Carlo simulations (histogram) and from the Planck 2015 maps (vertical pink lines) outside the UPB77 mask. Columns from left to right show different resolutions (N_{side} = 1024, 256, and 64), while rows show results for the four component separation methods. Unlike in Fig. 16, the variance distributions are not normalized here. 

In the text 
Fig. 18 Difference between the Npoint functions for the highpass filtered N_{side} = 64Planck 2015 CMB estimates and the corresponding means estimated from 1000 Monte Carlo simulations. The Stokes parameters Q_{r} and U_{r} were locally rotated so that the correlation functions are independent of coordinate frame. The first row shows results for the 2point function, from left to right, TQ_{r}, TU_{r}, Q_{r}Q_{r}, and U_{r}U_{r}. The second row shows results for the pseudocollapsed 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}, and the third row shows results for the equilateral 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r} and U_{r}U_{r}U_{r}. The red solid, orange tripledotdashed, green dashed and blue dotdashed lines correspond to the Commander, NILC, SEVEM, and SMICA maps, respectively. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using SMICA simulations. See Sect. 7.2 for the definition of the separation angle θ. 

In the text 
Fig. 19 Lensinginduced Bmode power spectra in the componentseparated polarization CMB maps. The solid line represents the bestfit cosmology from the Planck data release in 2015. Error bars were evaluated using a semianalytical approximation validated over the FFP8 simulations, as described in Planck Collaboration PIP116 (in prep.). 

In the text 
Fig. A.1 Commander processing masks for temperature (top) and polarization (bottom). For temperature, the different shades of grey correspond to different angular resolutions, ranging from 5′ (light grey) through 7.5′ (dark grey) to 40′ FWHM (black). For polarization, the same mask is used for both 10′ and 40′ FWHM resolution. 

In the text 
Fig. A.2 Multipole moment weights used for multiresolution hybridization in the Commander CMB map, as described by Eq. (A.1). 

In the text 
Fig. A.3 5° × 5° zoomin of the multiresolution contributions to the Commander hybrid CMB map from the 40′ (top left), (top right), and approximately 5′ (bottom left) solutions, centred on the south galactic pole. The hybrid map is shown in the bottom right panel. 

In the text 
Fig. A.4 Difference between the Commander CMB solutions derived using (1) only Planck data and (2) Planck+WMAP+408 MHz, both smoothed to a common resolution of 1° FWHM. The grey region indicates the common mask. 

In the text 
Fig. B.1 Needlet bands used in the analysis. The thick black line shows the normalization of the needlet bands, that is, the total filter applied to the original map after needlet decomposition and synthesis of the output map from needlet coefficients. 

In the text 
Fig. B.2 Fullsky average of needlet weights for different frequency channels and needlet bands. From top to bottom, the panels show results for temperature, Emodes, and Bmodes. 

In the text 
Fig. B.3 NILC masks for temperature (top) and polarization (bottom). 

In the text 
Fig. C.1 Weights used to combine the 143 and 217 GHz cleaned frequency channel maps into the final SEVEM temperature map. The weights do not sum to unity because they include the effect of deconvolving the beams of the frequency maps and convolving with the 5′ beam of the final map. 

In the text 
Fig. C.2 SEVEM masks in temperature (top) and polarization (bottom). 

In the text 
Fig. D.1 SMICA weights for temperature (top) and polarization (bottom). For readability, the values are shown for input maps in units of antenna temperature. The plot goes up to ℓ ≃ 3600, but the output maps are synthesized using all multipoles up to ℓ = 4000. For polarization, the thick solid lines show the contribution of input Emodes to the CMB Emodes and the thick dashed lines show the same for the Bmodes. The thin lines, all close to zero, show “crosscontributions” of input Emodes to the CMB Bmodes and vice versa. 

In the text 
Fig. D.2 SMICA masks in temperature (top) and polarization (bottom). 

In the text 
Fig. E.1 Pairwise differences between CMB temperature maps obtained from FFP8 simulations. Prior to differencing, the maps have been smoothed to 80′ FWHM and downgraded to N_{side} = 128. 

In the text 
Fig. E.2 Pairwise differences between CMB Q maps obtained from FFP8 simulations. Smoothing and degrading is as in Fig. E.1. 

In the text 
Fig. E.3 Pairwise differences between CMB U maps obtained from FFP8 simulations. Smoothing and degrading is as in Figs. E.1 and E.2. 

In the text 
Fig. E.4 Difference between the Npoint functions for the highpass filtered N_{side} = 64 FFP8 CMB estimates and the corresponding means estimated from 1000 MC simulations. The Stokes parameters Q_{r} and U_{r} were locally rotated so that the correlation functions are independent of coordinate frame. The first row shows results for the temperature map, from left to right, results for the 2point, pseudocollapsed 3point, equilateral 3point and connected rhombic 4point functions. We note that the lines lie on top of each other. The second row shows results for the 2point function, from left to right, TQ_{r}, TU_{r}, Q_{r}Q_{r}, and U_{r}U_{r}. The third row shows results for the pseudocollapsed 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}, and the fourth row shows results for the equilateral 3point function, from left to right, TQ_{r}Q_{r}, TU_{r}U_{r}, Q_{r}Q_{r}Q_{r}, and U_{r}U_{r}U_{r}. The black solid, red tripledotdashed, orange dashed, green dotdashed, and blue long dashed lines correspond to the true, Commander, NILC, SEVEM, and SMICA maps, respectively. The true CMB map was analysed with added noise corresponding to the SMICA component separation method. The shaded dark and light grey regions indicate the 68% and 95% confidence regions, respectively, estimated using SMICA simulations. See Sect. 7.2 for the definition of the separation angle θ. 

In the text 
Fig. E.5 Power spectra of the foregroundcleaned CMB maps from FFP8 simulations. Top: TT power spectra evaluated using the FFP8UT74 mask. Bottom: EE power spectra evaluated using the FFP8UP76 mask. Thick lines show the spectra of signal plus noise estimated from the halfmission halfsum maps; thin lines show the noise levels from halfmission halfdifference maps. The black line shows the input spectrum. 

In the text 
Fig. E.6 TT angular power spectra of residuals from the indicated FFP8 components in the Planck 2015 CMB maps, compared with the predicted signal from the best fit cosmology. “Other” is the sum of CO, freefree, thermal and kinetic SZ, spinning dust, and synchrotron emission. The horizontal axis is linear in ℓ^{ 0.5}. 

In the text 
Fig. E.7 EE angular power spectra of residuals from the indicated FFP8 components in the Planck 2015 CMB maps, compared with the predicted signal from the best fit cosmology. “Other” is the sum of CO, freefree, thermal and kinetic SZ, spinning dust, CIB fluctuations, and unresolved radio and infrared sources. The horizontal axis is linear in ℓ^{ 0.5}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.