Issue 
A&A
Volume 571, November 2014
Planck 2013 results



Article Number  A12  
Number of page(s)  31  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201321580  
Published online  29 October 2014 
Planck 2013 results. XII. Diffuse component separation
^{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}
Aalto University Metsähovi Radio Observatory and Dept of Radio
Science and Engineering, PO Box
13000, 00076
Aalto,
Finland
^{3}
African Institute for Mathematical Sciences, ,
6–8 Melrose Road, Muizenberg,
Cape Town
7945, South
Africa
^{4}
Agenzia Spaziale Italiana Science Data Center, via del Politecnico
snc, 00133
Roma,
Italy
^{5}
Agenzia Spaziale Italiana, Viale Liegi 26,
00198
Roma,
Italy
^{6}
Astrophysics Group, Cavendish Laboratory, University of
Cambridge, J J Thomson
Avenue, Cambridge
CB3 0HE,
UK
^{7}
Astrophysics & Cosmology Research Unit, School of
Mathematics, Statistics & Computer Science, University of
KwaZuluNatal, Westville Campus,
Private Bag X54001, 4000
Durban, South
Africa
^{8}
Atacama Large Millimeter/submillimeter Array, ALMA Santiago
Central Offices, Alonso de Cordova
3107, Vitacura, Casilla
763 0355, Santiago, Chile
^{9}
CITA, University of Toronto, 60 St. George St., Toronto, ON
M5S 3H8,
Canada
^{10}
CNR – ISTI, Area della Ricerca, via G. Moruzzi 1,
56124
Pisa,
Italy
^{11}
CNRS, IRAP, 9
Av. colonel Roche, BP
44346, 31028
Toulouse Cedex 4,
France
^{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), 28006
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 Física Fundamental, Facultad de Ciencias,
Universidad de Salamanca, 37008
Salamanca,
Spain
^{21}
Departamento de Física, Universidad de Oviedo, ,
Avda. Calvo Sotelo s/n,
33007
Oviedo,
Spain
^{22}
Departamento de Matemáticas, Estadística y Computación,
Universidad de Cantabria, Avda. de
los Castros s/n, 39005
Santander,
Spain
^{23}
Department of Astronomy and Astrophysics, University of
Toronto, 50 Saint George Street,
Toronto, Ontario,
Canada
^{24}
Department of Astrophysics/IMAPP, Radboud University
Nijmegen, PO Box
9010, 6500 GL
Nijmegen, The
Netherlands
^{25}
Department of Electrical Engineering and Computer Sciences,
University of California, Berkeley, California, USA
^{26}
Department of Physics & Astronomy, University of British
Columbia, 6224 Agricultural Road,
Vancouver, British
Columbia, Canada
^{27}
Department of Physics and Astronomy, Dana and David Dornsife
College of Letter, Arts and Sciences, University of Southern California, ,
Los Angeles, CA
90089,
USA
^{28}
Department of Physics and Astronomy, University College
London, London
WC1E 6BT,
UK
^{29}
Department of Physics, Florida State University, ,
Keen Physics Building, 77 Chieftan
Way, Tallahassee,
Florida,
USA
^{30}
Department of Physics, Gustaf Hällströmin katu 2a, University of
Helsinki, 00014
Helsinki,
Finland
^{31}
Department of Physics, Princeton University, ,
Princeton, New Jersey, USA
^{32}
Department of Physics, University of California, ,
One Shields Avenue,
Davis, California, USA
^{33}
Department of Physics, University of California, ,
Santa Barbara, California, USA
^{34}
Department of Physics, University of Illinois at
UrbanaChampaign, 1110 West Green
Street, Urbana,
Illinois,
USA
^{35}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli
Studi di Padova, via Marzolo
8, 35131
Padova,
Italy
^{36}
Dipartimento di Fisica e Astronomia, Università degliStudi di
Bologna, viale Berti Pichat
6/2, 40127
Bologna,
Italy
^{37}
Dipartimento di Fisica e Scienze della Terra, Università di
Ferrara, via Saragat
1, 44122
Ferrara,
Italy
^{38}
Dipartimento di Fisica, Università La Sapienza, ,
P.le A. Moro 2, 00185
Roma,
Italy
^{39}
Dipartimento di Fisica, Università degli Studi di
Milano, via Celoria,
16, 20133
Milano,
Italy
^{40}
Dipartimento di Fisica, Università degli Studi di
Trieste, via A. Valerio
2, 34127
Trieste,
Italy
^{41}
Dipartimento di Fisica, Università di Roma Tor
Vergata, via della Ricerca
Scientifica 1, 00133
Roma,
Italy
^{42}
Discovery Center, Niels Bohr Institute, ,
Blegdamsvej 17, 2100
Copenhagen,
Denmark
^{43}
Dpto. Astrofísica, Universidad de La Laguna (ULL), ,
38206, La Laguna, Tenerife, Spain
^{44}
European Southern Observatory, ESO Vitacura, Alonso de Cordova 3107, Vitacura, Casilla
19001
Santiago,
Chile
^{45}
European Space Agency, ESAC, Planck Science Office, Camino bajo
del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la
Cañada, 28692
Madrid,
Spain
^{46}
European Space Agency, ESTEC, Keplerlaan 1,
2201 AZ
Noordwijk, The
Netherlands
^{47}
Haverford College Astronomy Department, ,
370 Lancaster Avenue,
Haverford, Pennsylvania, USA
^{48}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2,
University of Helsinki, 00014
Helsinki,
Finland
^{49}
INAF – Osservatorio Astrofisico di Catania, via S. Sofia
78, 95123
Catania,
Italy
^{50}
INAF – Osservatorio Astronomico di Padova, Vicolo
dell’Osservatorio 5, 35122
Padova,
Italy
^{51}
INAF – Osservatorio Astronomico di Roma, via di Frascati
33, 00040
Monte Porzio Catone,
Italy
^{52}
INAF – Osservatorio Astronomico di Trieste, via G.B. Tiepolo
11, 34131
Trieste,
Italy
^{53}
INAF Istituto di Radioastronomia, via P. Gobetti
101, 40129
Bologna,
Italy
^{54}
INAF/IASF Bologna, via Gobetti 101, 40129
Bologna,
Italy
^{55}
INAF/IASF Milano, via E. Bassini 15, 20133
Milano,
Italy
^{56}
INFN, Sezione di Bologna, via Irnerio 46,
40126
Bologna,
Italy
^{57}
INFN, Sezione di Roma 1, Università di Roma
Sapienza, Piazzale Aldo Moro
2, 00185
Roma,
Italy
^{58}
INFN/National Institute for Nuclear Physics, ,
via Valerio 2, 34127
Trieste,
Italy
^{59}
IPAG: Institut de Planétologie et d’Astrophysique de Grenoble,
Université Joseph Fourier, Grenoble
1/CNRSINSU, UMR
5274, 38041
Grenoble,
France
^{60}
ISDC Data Centre for Astrophysics, University of
Geneva, ch. d’Ecogia
16, 1290
Versoix,
Switzerland
^{61}
IUCAA, Post Bag 4, Ganeshkhind, Pune University
Campus, 411 007
Pune,
India
^{62}
Imperial College London, Astrophysics group, Blackett
Laboratory, Prince Consort
Road, London,
SW7 2AZ,
UK
^{63}
Infrared Processing and Analysis Center, California Institute of
Technology, Pasadena,
CA
91125,
USA
^{64}
Institut Néel, CNRS, Université Joseph Fourier Grenoble
I, 25 rue des
Martyrs, 38042
Grenoble,
France
^{65}
Institut Universitaire de France, 103 Bd SaintMichel, 75005
Paris,
France
^{66}
Institut d’Astrophysique Spatiale, CNRS (UMR8617) Université
ParisSud 11, Bâtiment
121, 91405
Orsay,
France
^{67}
Institut d’Astrophysique de Paris, CNRS (UMR7095), ,
98bis Boulevard Arago,
75014
Paris,
France
^{68}
Institute for Space Sciences, 077125
BucharestMagurale,
Romania
^{69}
Institute of Astronomy and Astrophysics, Academia
Sinica, 0617
Taipei,
Taiwan
^{70}
Institute of Astronomy, University of Cambridge, ,
Madingley Road, Cambridge
CB3 0HA,
UK
^{71}
Institute of Theoretical Astrophysics, University of
Oslo, Blindern,
0315
Oslo,
Norway
^{72}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38205,
La Laguna, Tenerife,
Spain
^{73}
Instituto de Física de Cantabria (CSICUniversidad de
Cantabria), Avda. de los Castros
s/n, 39005
Santander,
Spain
^{74}
Istituto di Fisica del Plasma, CNRENEAEURATOM Association, via
R. Cozzi 53, 20125
Milano,
Italy
^{75}
Jet Propulsion Laboratory, California Institute of
Technology, 4800 Oak Grove
Drive, Pasadena,
California,
USA
^{76}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School
of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13
9PL, UK
^{77}
Kavli Institute for Cosmology Cambridge, ,
Madingley Road, Cambridge, CB3 0HA, UK
^{78}
Kavli Institute for Theoretical Physics, University of
California, Santa Barbara Kohn
Hall, Santa
Barbara, CA
93106,
USA
^{79}
LAL, Université ParisSud, CNRS/IN2P3, ,
91898
Orsay,
France
^{80}
LERMA, CNRS, Observatoire de Paris, 61 Avenue de
l’Observatoire, Paris, France
^{81}
Laboratoire AIM, IRFU/Serviced’Astrophysique – CEA/DSM – CNRS –
Université Paris Diderot, Bât. 709, CEASaclay, 91191
GifsurYvette Cedex,
France
^{82}
Laboratoire Traitement et Communication de l’Information, CNRS
(UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634
Paris Cedex 13,
France
^{83}
Laboratoire de Physique Subatomique et de Cosmologie, Université
Joseph Fourier Grenoble I, CNRS/IN2P3, Institut National Polytechnique de
Grenoble, 53 rue des
Martyrs, 38026
Grenoble Cedex,
France
^{84}
Laboratoire de Physique Théorique, Université ParisSud 11
& CNRS, Bâtiment
210, 91405
Orsay,
France
^{85}
Lawrence Berkeley National Laboratory, ,
Berkeley, California, USA
^{86}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741
Garching,
Germany
^{87}
McGill Physics, Ernest Rutherford Physics Building, McGill
University, 3600 rue University, Montréal, QC,
H3A 2T8,
Canada
^{88}
MilliLab, VTT Technical Research Centre of Finland, Tietotie
3, 02044
Espoo,
Finland
^{89}
National University of Ireland, Department of Experimental
Physics, Maynooth,
Co. Kildare,
Ireland
^{90}
Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen,
Denmark
^{91}
Observational Cosmology, Mail Stop 36717, California Institute of
Technology, Pasadena,
CA, 91125, USA
^{92}
Optical Science Laboratory, University College
London, Gower
Street, London,
UK
^{93}
SBITPLPPC, EPFL, 1015, Lausanne,
Switzerland
^{94}
SISSA, Astrophysics Sector, via Bonomea 265,
34136
Trieste,
Italy
^{95}
School of Physics and Astronomy, Cardiff University, ,
Queens Buildings, The Parade,
Cardiff, CF24 3AA, UK
^{96}
School of Physics and Astronomy, University of
Nottingham, Nottingham
NG7 2RD,
UK
^{97}
Space Research Institute (IKI), Russian Academy of
Sciences, Profsoyuznaya Str,
84/32, 117997
Moscow,
Russia
^{98}
Space Sciences Laboratory, University of California, ,
Berkeley, California, USA
^{99}
Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo
Mall, Stanford,
California,
USA
^{100}
SubDepartment of Astrophysics, University of
Oxford, Keble Road,
Oxford
OX1 3RH,
UK
^{101}
Theory Division, PHTH, CERN, 1211, Geneva 23, Switzerland
^{102}
UPMC Univ Paris 06, UMR7095, 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 Infrared Astronomy, 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}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478
Warszawa,
Poland
Received: 26 March 2013
Accepted: 18 December 2013
Planck has produced detailed allsky observations over nine frequency bands between 30 and 857 GHz. These observations allow robust reconstruction of the primordial cosmic microwave background (CMB) temperature fluctuations over nearly the full sky, as well as new constraints on Galactic foregrounds, including thermal dust and line emission from molecular carbon monoxide (CO). This paper describes the component separation framework adopted by Planck for many cosmological analyses, including CMB power spectrum determination and likelihood construction on large angular scales, studies of primordial nonGaussianity and statistical isotropy, the integrated SachsWolfe effect, gravitational lensing, and searches for topological defects. We test four foregroundcleaned CMB maps derived using qualitatively different component separation algorithms. The quality of our reconstructions is evaluated through detailed simulations and internal comparisons, and shown through various tests to be internally consistent and robust for CMB power spectrum and cosmological parameter estimation up to ℓ = 2000. The parameter constraints on ΛCDM cosmologies derived from these maps are consistent with those presented in the crossspectrum based Planck likelihood analysis. We choose two of the CMB maps for specific scientific goals. We also present maps and frequency spectra of the Galactic lowfrequency, CO, and thermal dust emission. The component maps are found to provide a faithful representation of the sky, as evaluated by simulations, with the largest bias seen in the CO component at 3%. For the lowfrequency component, the spectral index varies widely over the sky, ranging from about β = −4 to − 2. Considering both morphology and prior knowledge of the low frequencycomponents, the index map allows us to associate a steep spectral index (β< −3.2) with strong anomalous microwave emission, corresponding to a spinning dust spectrum peaking below 20 GHz, a flat index of β> −2.3 with strong freefree emission, and intermediate values with synchrotron emission.
Key words: cosmic background radiation
© ESO, 2014
1. Introduction
This paper, one of a set associated with the 2013 release of data from the Planck^{1} mission (Planck Collaboration I 2014), describes the component separation techniques applied to the Planck data to produce maps of the cosmic microwave background (CMB) temperature anisotropies (see Fig. 1) and of diffuse foregrounds.
Fig. 1
Foregroundcleaned CMB maps derived by CommanderRuler, NILC, SEVEM and SMICA. Note that the SMICA map has been filled in smoothly inside a 3% Galactic mask. 
The sky at millimetre and submillimetre wavelengths contains a wealth of cosmological and astrophysical information. Accessing it is an inversion process, known as component separation, to extract the sources of emission contributing to a set of maps observed at different frequencies. Planck gives us a powerful data set to unlock new information in this manner by observing the entire sky from 30 to 857 GHz in nine frequency bands at higher angular resolution and sensitivity than its predecessors. Accurate and detailed component separation is a central objective of the mission.
We divide the foregrounds into two distinct categories: diffuse emission from the Galaxy and compact sources. The Galactic foregrounds are the principal source of contamination of the CMB on large angular scales, with fluctuation power decreasing roughly as a power law towards higher multipoles (Bennett et al. 2003). They are dominated by synchrotron, freefree and anomalous microwave emission (AME, ascribed to spinning dust grains) at frequencies below 70 GHz, and by rotational line emission from carbon monoxide (CO) molecules and thermal dust emission at frequencies above 100 GHz. Extragalactic foregrounds, on the other hand, dominate the smallscale contamination of the CMB. They arise from discrete, individually detectable compact sources and the collective emission from unresolved radio and infrared (IR) sources, and also from the SunyaevZeldovich (SZ) effect in galaxy clusters (Planck Collaboration XVIII 2011; Planck Collaboration XXVIII 2014; Planck Collaboration XXIX 2014).
In the Planck analyses, these foregrounds are dealt with in a variety of ways. At the power spectrum and likelihood level, the extragalactic foregrounds are modelled with parameterized power spectra, appropriate to their statistical isotropy, over regions restricted to low Galactic emission (Planck Collaboration XV 2014). Component separation as described in the present paper aims at removing Galactic emission to produce CMB maps covering the largest possible sky area for studies of the largescale properties and higherorder statistics of the CMB. In addition, this component separation provides a reconstruction of the diffuse emission from our Galaxy. Detailed studies of specific extragalactic foregrounds, such as the cosmic infrared background (CIB; Planck Collaboration XVIII 2014) and the diffuse SZ signal (Planck Collaboration XXI 2014), employ methods tailored to their particular needs.
Building on previous work (Leach et al. 2008), we approach CMB extraction with a philosophy designed to ensure robustness by applying four distinct algorithms based on two different methodologies. The first avoids any assumptions concerning the foregrounds and relies solely on a minimum variance criterion for the data component possessing a blackbody spectrum (i.e., the CMB), while the second methodology relies on parametric modelling of the foregrounds in either real or harmonic space. We evaluate the performance of these component separation algorithms through detailed simulations, and we examine the robustness of the recovered CMB maps by comparing them, their power spectra, and their resulting cosmological constraints. As a diagnostic, we also briefly examine their higherorder statistics.
The CMB results presented in this work serve a number of applications. We use the realspace modelling to produce a clean CMB map and power spectra on large angular scales, where diffuse Galactic emission is the main contaminant, to construct the likelihood function at low multipoles; this is then combined with the high multipole likelihood function that models extragalactic foregrounds with power spectra (Planck Collaboration XV 2014). The high resolution CMB maps are used as a check on primary cosmological constraints (see below), for lensing studies (Planck Collaboration XVII 2014), studies of the integrated SachsWolfe ISW effect (Planck Collaboration XIX 2014), of the isotropy of the CMB (Planck Collaboration XXIII 2014), of nonGaussian statistics (Planck Collaboration XXIV 2014), in searches for topological defects (Planck Collaboration XXV 2014), and for examination of the geometry and topology of the Universe (Planck Collaboration XXVI 2014).
In addition, we present maps of diffuse Galactic emission divided into low and highfrequency components, as well as a molecular CO component. We judge the adequacy of this reconstruction through simulations and by comparison with known properties of the diffuse Galactic foregrounds.
The paper is organized as follows. In Sect. 2 we discuss the expected sources of sky emission over the Planck frequency interval and how they are modelled. Then in Sect. 3 we detail the overall approach and introduce the four component separation methods. In Sect. 4 we present the Planck data set and preprocessing procedure, and we describe our simulations. This is followed by a presentation of the derived CMB maps and their characterization in Sect. 5. Section 6 is dedicated to power spectra and cosmological parameter constraints obtained from these maps, and Sect. 7 to studies of higherorder statistics. Section 8 presents a reconstruction of the diffuse Galactic foregrounds, and Sect. 9 concludes. We relegate details of the algorithms to appendices.
2. The sky at Planck frequencies
The properties of Galactic emission vary significantly across the Planck frequency range from 30 to 857 GHz. At frequencies below 70 GHz, the dominant radiation processes are: synchrotron emission from cosmic ray electrons interacting with the Galactic magnetic field (e.g., Haslam et al. 1982; Reich & Reich 1988; Broadbent et al. 1989; Davies et al. 1996; Platania et al. 2003; Bennett et al. 2003; Gold et al. 2011); thermal Bremsstrahlung (or freefree emission) from electronelectron and electronion scattering (e.g., Banday et al. 2003; Dickinson et al. 2003; Davies et al. 2006; Ghosh et al. 2012; Planck Collaboration Int. XII 2013; Planck Collaboration XX 2011); and AME from dust grains (Kogut 1996; Leitch et al. 1997; Banday et al. 2003; Lagache 2003; de OliveiraCosta et al. 2004; Finkbeiner et al. 2004; Davies et al. 2006; Bonaldi et al. 2007; Dobler & Finkbeiner 2008; MivilleDeschênes et al. 2008; Ysard et al. 2010; Gold et al. 2011; Planck Collaboration XX 2011), possibly due to their rotational emission (Draine & Lazarian 1998; AliHaïmoud et al. 2009; Ysard & Verstraete 2010; Hoang & Lazarian 2012). Over the frequency range covered by Planck, both synchrotron and freefree spectra are well approximated by power laws in antenna temperature, T_{B} ∝ ν^{β}, with the synchrotron index, β_{synch}, ranging from − 3.2 to − 2.8 (Davies et al. 1996) and the freefree index, β_{ff}, lying between − 2.2 and − 2.1. Less is known about the AME spectrum, but spinning dust models with a spectrum peaking at frequencies below 20 GHz (in antenna temperature units) adequately describe current observations^{2}. Above the peak, the spectrum appears consistent with a powerlaw (e.g., Banday et al. 2003; Davies et al. 2006; Dobler & Finkbeiner 2008; Ghosh et al. 2012). In addition to these three, the existence of a fourth lowfrequency foreground component, known as the “Galactic haze”, has been claimed, possibly due to a hardspectrum synchrotron population near the Galactic centre (e.g., Finkbeiner 2004; Dobler & Finkbeiner 2008; Pietrobon et al. 2012; Planck Collaboration Int. IX 2013).
At frequencies higher than 100 GHz, thermal dust emission dominates over most of the sky and is commonly described by a modified blackbody spectrum with powerlaw emissivity, ϵ_{ν} ∝ ν^{βd}, and temperature, T_{d}. Both the temperature and spectral index, β_{d}, vary spatially. Prior to Planck, the bestfitting single component dust model had a temperature T_{d} ≈ 18 K and spectral index β_{d} ≈ 1.7 (Finkbeiner et al. 1999; Bennett et al. 2003; Gold et al. 2011), although there is evidence of flattening of the spectral index from around 1.8 in the farinfrared to 1.55 in the microwave region (Finkbeiner et al. 1999; Planck Collaboration Int. XIV 2014), the interpretation of which is still under study.
In addition to these diffuse Galactic components, extragalactic emission contributes at Planck frequencies. In particular, a large number of radio and farinfrared (FIR; Planck Collaboration XIII 2011) galaxies, clusters of galaxies and the CIB (Planck Collaboration XVIII 2011) produce a statistically isotropic foreground, with frequency spectra well approximated by models similar to those applicable to the Galactic foregrounds (modified blackbody spectra, power laws, etc.). Except for a frequencydependent absolute offset, which may be removed as part of the overall offset removal procedure, these extragalactic components are therefore typically absorbed by either the lowfrequency or thermal dust components during component separation. No special treatment is given here to extragalactic foregrounds, beyond the masking of bright objects. Dedicated scientific analyses of these sources are described in detail in Planck Collaboration XVIII (2011), Planck Collaboration XXVIII (2014), and Planck Collaboration XXIX (2014). In the Planck likelihood, extragalactic sources are modelled in terms of power spectrum templates at high ℓ (Planck Collaboration XV 2014).
Other relevant sources include emission from molecular clouds, supernova remnants, and compact H ii regions inside our own Galaxy, as well as the thermal and kinetic SZ effects, due to inverse Compton scattering of CMB photons off free electrons in ionized media. Planck provides new and important information on all these processes, as described both in the following and in the companion papers Planck Collaboration XIII (2014), Planck Collaboration XXI (2014), and Planck Collaboration XI (2014). In particular, Planck’s frequency range, angular resolution and sensitivity make it a powerful probe of thermal dust, resulting in new and tight constraints on dust temperature and emissivity. The same frequencies also allow extraction of the first ever fullsky maps of the emission resulting from the CO J = 1 → 0, J = 2 → 1 and J = 3 → 2 rotational transitions at 115, 230 and 345 GHz, respectively (Planck Collaboration XIII 2014).
The focus of this paper is to reconstruct the CMB anisotropies over a large sky fraction, exploiting only the Planck frequency bands. We also present a detailed reconstruction of the thermal dust emission at high frequencies, as well as CO emission lines. At low frequencies and over the region used for CMB analysis, the total foreground contribution is well approximated by a single power law (see Sect. 8). We therefore model the sum of all lowfrequency foregrounds by a power law with spatially varying spectral index whose numerical value in any pixel results from the influence of the dominant foreground component at that location. The full analysis of diffuse foregrounds, using ancillary data to resolve the individual components at low frequencies, will be presented in a forthcoming publication.
Overview and comparison of component separation algorithms.
3. Approach to component separation
The rich content of the Planck data encourages application of several component separation techniques. We consider four, as summarized in Table 1, which we classify according to one of two different general methodologies. The first makes minimal assumptions concerning the foregrounds and seeks only to minimize the variance of the CMB, i.e., the sky component possessing a blackbody spectrum. We implement this approach with a needlet (wavelet on the sphere) version of the internal linear combination (ILC) algorithm (NILC; Delabrouille et al. 2009), and also with a templatebased method to remove foreground contamination from the CMBdominant bands. These foreground templates are constructed from the lowest and highest frequency channels (FernándezCobos et al. 2012, Spectral Estimation Via Expectation Maximization, SEVEM).
The second methodology uses parametric modelling of the foregrounds. In our real space implementation, we explore model parameters through Bayesian parameter estimation techniques, fitting a parametric signal model per pixel (Commander; Eriksen et al. 2006, 2008); a similar implementation is presented by Stompor et al. (2009). To estimate spectral indices robustly in pixel space, this procedure requires identical angular resolution across all frequencies included in the analysis, and is therefore limited in resolution by the 30 GHz LFI channel. However, this is sufficient to generate the lowresolution CMB map and power spectrum samples required for the low multipole part of the Planck likelihood function for cosmological parameters (Planck Collaboration XV 2014). To produce full resolution maps, we use the resulting lowresolution spectral parameter samples to solve for the component amplitudes, in an extension to the method known as Ruler (we refer to the combined method as CommanderRuler, or CR). In our fourth technique, we implement a CMBoriented parametric approach that fits the amplitude and spectral parameters of CMB and foregrounds in the harmonic domain (Spectral Matching Independent Component Analysis, SMICA; Cardoso et al. 2008).
Details of each algorithm are given in the appendices. We now turn to their application to the data and evaluate their performance using simulations.
4. Data, simulations and masks
We use the data set from the first 15.5 months of Planck observations, corresponding to 2.6 sky surveys, from both the Low Frequency Instrument (LFI) and High Frequency Instrument (HFI). The primary inputs for component separation are the frequency channel maps, including halfring maps, bandpasses, and beam characteristics; a full description of these products is given in Planck Collaboration II (2014) and Planck Collaboration VI (2014). No special corrections are made for zodiacal light emission (ZLE; Planck Collaboration VI 2014) in the analyses presented here. The ZLE is not stationary on the sky, since it depends on Planck’s position and scanning strategy. Therefore the frequency maps contain a projected version of the emission averaged over the nominal mission. Despite this, a series of exploratory analyses showed that our algorithms naturally correct for this component within their existing model space. It was also found that larger CMB residuals were induced when applying a correction based on a ZLE model than when applying no correction, most likely due to uncertainties in the model itself.
To evaluate and validate our algorithms, we analyse a large suite of realistic simulations, the socalled full focal plane (FFP) simulations, based on detailed models of the instrument and sky. The version used for this data release is denoted FFP6, and is described in Planck Collaboration (2013). The simulation procedure generates time streams for each detector, incorporating the satellite pointing, the individual detector beams, bandpasses, noise properties, and data flags, and then produces simulated frequency channel maps through the mapmaking process. For the input sky, we use the Planck Sky Model (PSM), which includes the CMB, diffuse Galactic emission (synchrotron, freefree, thermal dust, AME, and molecular CO lines), and compact sources (thermal and kinetic SZ effects, radio sources, infrared sources, the CIB, and ultracompact H ii regions). The prelaunch version of the PSM is described by Delabrouille et al. (2013), and has been modified for the present work as described in Planck Collaboration (2013). Each FFP data set consists of three parts: the simulated observations, Monte Carlo realizations of the CMB, and Monte Carlo realizations of the instrumental noise.
For both the data and the simulations, we reconstruct the CMB and foregrounds from the full frequency channel maps and the corresponding halfring maps, which are made from the data in the first half or second half of each stable pointing period. The halfring maps can be used to obtain an estimate of the noise in each channel by taking half of the difference between the two maps, thereby normalizing the noise level to that of the full map. This is referred to as the halfring halfdifference (HRHD) map. The signals fixed to the sky will be cancelled leaving only the noise contribution. The HRHD map can be treated as a realization of the same underlying noise processes and it can be used to estimate the power spectrum, and other properties, of the noise. If there are noise correlations between the halfring maps, then the estimates of the noise properties thus obtained can be biased. This is the case for HFI channels; the cosmic ray glitch removal (Planck Collaboration VI 2014; Planck Collaboration X 2014) induces correlations that lead to the noise power spectrum being underestimated by a few percent at high ℓ when using the HRHD maps.
Prior to processing the data through each component separation pipeline, we define masks for the point sources and bright Galactic regions. Point source masking is based on the source catalogues obtained by filtering the input sky maps with the Mexican Hat wavelet 2 (MHW2) filter and applying a 4σ threshold for the LFI bands and a 5σ threshold for the HFI bands (Planck Collaboration XVIII 2011; Planck Collaboration XXVIII 2014). The mask radius of each source is different for the LFI and HFI. Due to the large beam size of LFI channels, we define a variable masking radius for each source according to its signaltonoise ratio (S/N) as , where r is the radius, A is the S/N, and m is the maximum amplitude (given in units of the background noise level) allowed for the tail of unmasked point sources; we set m = 0.1, which is a compromise between masking the source tails and minimizing the number of masked pixels. For HFI, the mask radius around each source is 1.27 × FWHM, using the average FWHM obtained from the effective beams.
A basic set of Galactic masks is defined as follows. We subtract a CMB estimate from the 30 and 353 GHz maps, mask point sources, and smooth the resulting maps by a Gaussian with FWHM of 5°. We then threshold and combine them, generating a series of masks with different amounts of available sky. The resulting combined Galactic (CG) masks, shown in Fig. 2, correspond to sky fractions of 20, 40, 60, 70, 75, 80, 90, 97, and 99%, and are named CG20, etc.
Fig. 2
Combined Galactic (CG) emission masks for the Planck data, corresponding to sky fractions of 20, 40, 60, 70, 75, 80, 90, 97, and 99%. The masks are named CG20, etc. 
5. CMB Maps
We begin the discussion of our results by presenting the foregroundcleaned CMB maps. These maps are shown in Fig. 1 for each of the four component separation algorithms. Already from this figure it is clear that the wide frequency coverage and high angular resolution of Planck allow a faithful reconstruction of the CMB field over most of the sky. The fluctuations appear visually consistent with the theoretical expectation of a Gaussian and isotropic signal everywhere except inside a small band very close to the Galactic plane^{3}
Each CMB map is accompanied by its own confidence mask outside which the corresponding solution is considered statistically robust, shown in Fig. 3; for a definition of each mask, see Appendices A–D. Accepted sky fractions are 75, 93, 76, and 89%, respectively, for CommanderRuler, NILC, SEVEM, and SMICA. These masks are denoted CSCR75, CSNILC93, CSSEVEM76, CSSMICA89, respectively. The union of the confidence masks accepts 73% of the sky and is denoted U73. It is adopted as the default mask for evaluation purposes in this paper.
In addition to the CMB maps from the full data set, the halfring frequency maps have been processed by each algorithm to provide halfring CMB maps. They are used to provide estimates of the instrumental noise contribution to the foregroundcleaned maps in the power spectrum analysis (see Sect. 6). The algorithms were also used to process Monte Carlo simulations: 1000 realizations of the CMB and 1000 realizations of noise. They are not used in the analyses presented in this paper, but are used by Planck Collaboration XXIII (2014) and Planck Collaboration XXIV (2014).
Fig. 3
Summary of component separation (CS) confidence masks. Each pixel is encoded in terms of a sum in which CommanderRuler equals 1 (light blue), NILC equals 2 (dark red), SEVEM equals 4 (yellow), and SMICA equals 8 (light red). The masks are named CSCR75, CSNILC93, CSSEVEM76, and CSSMICA89, respectively, reflecting their accepted sky fraction. The union mask (U73), used for evaluation purposes in this paper, removes all coloured pixels. 
Fig. 4
Beam transfer functions of the four foregroundcleaned CMB maps. 
Fig. 5
Standard deviation between the four foregroundcleaned CMB maps. All maps have been downgraded to a HEALPix resolution of N_{side} = 128. The differences are typically less than 5 μK at high Galactic latitudes, demonstrating that the maps are consistent over a large part of the sky. 
The beam transfer functions of the foregroundcleaned CMB maps have been estimated for each algorithm, as shown in Fig. 4. The angular resolution of the NILC, SEVEM, and SMICA maps corresponds to a Gaussian beam with FWHM of 5′. The difference between SEVEM and NILC/SMICA is due to their different treatment of the HEALPix^{4} pixel window function (Górski et al. 2005). The deviation of NILC beam from a Gaussian shape at ℓ > 2800 is caused by the last needlet window (see Appendix B). CommanderRuler has a larger beam, because it is defined explicitly as a weighted average of frequency maps in pixel space. Its resolution is equivalent to a Gaussian beam with FWHM of approximately . The beam transfer functions have been computed assuming the bestfit beam transfer function for each frequency channel, and the uncertainties in the latter have not been propagated to these estimates.
Fig. 6
Pairwise differences between foregroundcleaned CMB maps. All maps have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale differences. The linelike discontinuities in the differences involving SEVEM is due to the two different regions used in this algorithm to clean the sky (see Appendix C for details). 
Fig. 7
CMB residual maps from the FFP6 simulation. A monopole determined at high Galactic latitude has been subtracted from the maps, and they have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale features. The residuals presented here provide a conservative estimate of those expected in the data (see text for details). 
In Fig. 5 we show the standard deviation per pixel among the four foregroundcleaned CMB maps downgraded to N_{side} = 128, and in Fig. 6 we show all pairwise difference maps. Typical differences at high Galactic latitudes are smaller than 5 μK. Considering the difference maps in more detail, it is clear that the CommanderRuler map is the most different from the other three, whereas NILC and SMICA are the most similar. This is not completely unexpected, because while CommanderRuler uses only frequencies between 30 and 353 GHz in its solution, the other three codes additionally include the dustdominated 545 and 857 GHz maps.
This difference in data selection may explain some of the coherent structures seen in Fig. 6. In particular, the most striking largescale feature in the difference maps involving CommanderRuler is a large negative band roughly following the ecliptic plane. This is where the ZLE (Planck Collaboration VI 2014) is brightest. Since the ZLE is also stronger at high frequencies, having a spectrum close to that of thermal dust, it is possible that this pattern may be an imprint of residual ZLE either in the CommanderRuler map, or in all of the other three maps. Both cases are plausible. The CommanderRuler solution may not have enough highfrequency information to distinguish between ZLE and normal thermal dust emission, and, by assuming a thermal dust spectrum for the entire highfrequency signal at 353 GHz, oversubtracts the ZLE at lower frequencies. It is also possible that the other three CMB solutions have positive ZLE residuals from extrapolating the highfrequency signal model from 857 GHz to the CMB frequencies. Without an accurate and detailed ZLE model, it is difficult to distinguish between these two possibilities. It is of course also possible that the true explanation is in fact unrelated to ZLE, and the correlation with the ecliptic plane is accidental. In either case, it is clear that the residuals are small in amplitude, with peaktopeak values typically smaller than 10 μK, of which by far the most is contained in a quadrupole aligned with the ecliptic. This provides additional evidence that residual ZLE is not important for the CMB power spectrum and cosmological parameter estimation, although some care is warranted when using these maps to study the statistics of the very largest angular scales (e.g., Planck Collaboration XXIII 2014); checking consistency among all four maps for a given application alleviates much of this concern.
We end this section by showing in Fig. 7 a set of residual maps derived by analysing the FFP6 simulation with exactly the same analysis approaches as applied to the data. It is evident that SMICA produces the map with lowest level of residuals. Considering the morphology in each case, we see that the main contaminant for CommanderRuler is undersubtracted freefree emission, while for both NILC and SEVEM it is oversubtracted thermal dust emission, and for SMICA it is undersubtracted thermal dust emission. However, at high latitudes and outside the confidence masks, the residuals are generally below a few μK in amplitude. It is also worth noting that each algorithm has been optimized (in terms of model definition, localization parameters, etc.) for the data, and the same configuration was subsequently used for the FFP6 simulations without further tuning. The simulations presented here therefore provide a conservative estimate of the residuals in the data. This is also reflected in the fact that the differences between CMB reconstructions for the FFP6 simulations are larger than those found in the data. See Appendix E for further details.
6. Power spectrum and cosmological parameters
In this section we evaluate the foregroundcleaned maps in terms of CMB power spectra and cosmological parameters. Our purpose in doing this is to show that the maps are consistent with the highℓ likelihood obtained from the crossspectrum analysis of detector set and frequency maps in Planck Collaboration XV (2014), and with the cosmological parameters derived from them in Planck Collaboration XVI (2014). This also establishes the consistency between Planck’s cosmological constraints and studies of the largescale structure and higher order statistics of the CMB.
6.1. Power spectra
Figure 8 shows the power spectra of the foregroundcleaned CMB maps and the corresponding HRHD maps, evaluated using the U73 mask with a 30′ cosine apodization. The spectra have been corrected for the effect of the mask and the beam transfer function of each algorithm has been deconvolved. The spectra of the HRHD maps give an estimate of the instrumental noise contribution to the power spectrum of the cleaned map. The correlations between the HFI halfring frequency maps are inherited by the halfring CMB maps that use them as input. At small angular scales, the CMB solution comes almost entirely from data in the HFI channels, and therefore the spectrum of the CMB HRHD maps is also biased low.
Fig. 8
Angular power spectra of the foregroundcleaned CMB maps and halfring halfdifference (HRHD) maps. The spectra have been evaluated using the U73 mask apodized with a 30′ cosine function. 
Fig. 9
Angular power spectra of FFP6 simulated components evaluated over the common mask (top) and the common point source mask (bottom), both apodized with a 30′ cosine function. Three components are shown: the CMB (dashed line); noise (dotdashed line); and the sum of all foregrounds (solid line). A nonlinear scale is used on the horizontal axis to show all the features of the spectra. 
At small angular scales, the effective noise levels of NILC, SEVEM, and SMICA are very similar, and lower than that of CommanderRuler. The last has larger noise because it operates entirely in pixel space and therefore applies the same weights to all multipoles. It cannot take advantage of the changing signaltonoise ratio of the frequency channels with angular scale.
We can estimate the contribution of residual foregrounds to the foregroundcleaned CMB maps by making use of the FFP6 simulations. In addition to processing the simulated frequency maps, the maps of the individual input sky components were processed by the algorithms after fixing their parameters or weights to the values obtained from the “observed” maps. Figure 9 shows the power spectra of the simulated FFP6 components, in this case CMB, noise and the sum of the foreground components. The top panel shows the spectra computed using the union mask derived from the simulation with a 30′ cosine apodization. The total foreground contribution becomes comparable to the CMB signal at ℓ ≈ 2000. The bottom panel shows the same computed with an apodized point source mask applied to the maps (i.e., no diffuse masking, although this mask does removes a large part of the Galactic plane). The residual foreground contribution is larger at all angular scales, but still it only becomes comparable to the CMB signal at ℓ ≈ 1800 in the worst case. For both masks, SMICA has the smallest residual foreground contamination at large angular scales, which is also demonstrated in Fig. 7. A more detailed examination of the contribution of the individual foreground components to the power spectrum is in Appendix E.
6.2. Likelihood and cosmological parameters
We estimate the binned power spectra with XFaster (Rocha et al. 2010, 2011) and determine cosmological parameter constraints using a correlated Gaussian likelihood. Parameter constraints are derived using a MetropolisHastings Markov Chain Monte Carlo sampler. To speed up this process, we additionally use PICO (Parameters for the Impatient COsmologist; Fendt & Wandelt 2008), a tool which interpolates the CMB power spectra and matter power spectra as a function of cosmological parameters.
6.2.1. Model and methods
We compute the power spectrum for each foregroundcleaned map over the multipole range 2 ≤ ℓ ≤ 2500, while parameter constraints are derived using only 70 ≤ ℓ ≤ 2000; as shown in Appendix E through simulations, modelling errors become nonnegligible between ℓ = 2000 and 2500. For parameter estimation, we adopt a standard sixparameter ΛCDM model, and impose an informative Gaussian prior of τ = 0.0851 ± 0.014, since polarization data are not included in this analysis.
While the foregroundcleaned maps should have minimal contamination from diffuse Galactic emission, they do contain significant contamination from unresolved extragalactic sources. These contributions are most easily modelled in terms of residual power spectra, therefore we marginalize over the corresponding parameters at the power spectrum level. To the six ΛCDM parameters, describing the standard cosmology, we add two foreground parameters, A_{ps}, the amplitude of a Poisson component (and hence constant, C_{ℓ} = A_{ps}), and A_{cl}, the amplitude of a clustered component with shape D_{ℓ} = ℓ(ℓ + 1)C_{ℓ}/ 2π ∝ ℓ^{ 0.8}. Both are expressed in terms of D_{ℓ} at ℓ = 3000 in units of μK^{2}.
The power spectrum calculation is based on the halfring halfsum (HRHS) and HRHD CMB maps (see Sect. 5); the latter is used to estimate the noise bias in the power spectra extracted from the HRHS maps. From these, we calculate the pseudospectra, and (Hivon et al. 2002), respectively, after applying the U73 mask. These are used as inputs to XFaster together with the beam transfer functions provided by each method (see Fig. 4).
To avoid aliasing of power from large to small scales, which would add an offset between the signalplusnoise and noise pseudospectra at high ℓ, we use the apodized version of the U73 mask. The known mismatch in the noise level between the spectra due to the correlation between the halfring maps is not explicitly corrected. It is left to be absorbed into the two foreground parameters.
Using the pseudospectra and XFaster, we then reconstruct an estimate of the power spectrum of each foregroundcleaned HRHS map, removing the noise bias as estimated from the corresponding HRHD map. To this end we apply an iterative scheme starting from a flat spectrum model. The result is a binned power spectrum and the associated Fisher matrix, which are then used to construct the likelihood, approximated here by a correlated Gaussian distribution.
To study consistency in the lowℓ range, we fit a twoparameter q–n (amplitudetilt) model relative to the Planck bestfit ΛCDM model on the form, , using a pixelspace likelihood for maps smoothed to 6° FWHM; see Planck Collaboration XV (2014) for further algorithmic details.
Fig. 10
Estimates of the CMB power spectra from the foregroundcleaned maps, computed by XFaster. The solid lines show the spectra after subtracting the bestfit model of residual foregrounds. The vertical dotted line shows the maximum multipole (ℓ = 2000) used in the likelihood for fitting the foreground model and cosmological parameters (see Sect. 6.2.2 for further details). The dashed lines show the spectra before residual foreground subtraction. 
Fig. 11
Comparison of cosmological and foreground parameter values estimated from the foregroundcleaned CMB maps for ℓ_{max} = 2000 (in red) and those obtained with CamSpec and Plik likelihoods (in blue). The values of the foreground parameters are not shown for CamSpec and Plik, since they use a different foreground model. 
6.2.2. Results
We perform the power spectrum and parameter estimation analysis for both the data and the FFP6 simulations described in Sect. 4. The results for the latter are given in Appendix E.
Figure 10 shows estimates of the angular power spectrum for each foregroundcleaned map, with the uncertainties given by the Fisher matrix. The parameter summary given in Fig. 11 shows the parameter constraints derived using multipoles between ℓ = 70 and 2000, and compares these to results obtained with the CamSpec and Plik likelihoods (Planck Collaboration XV 2014).
Differences in the power spectra at high ℓ are mostly absorbed by the twoparameter foreground model, rendering consistent cosmological parameters. For example, the highℓ power excess seen in the CommanderRuler map is wellfitted in terms of residual point sources, which makes intuitive sense, considering the lower angular resolution of this map (see Sect. 5). However, the ΛCDM parameter uncertainties derived from the four codes are very consistent. This indicates that most of the cosmological information content above ℓ ≥ 1500 is degenerate with the extragalactic foreground model, and a more sophisticated foreground treatment is required in order to recover significant cosmological parameter constraints from these scales. Beyond this, deviations among cosmological parameters are small and within 1σ for all methods and most of the parameters. Further, the parameters derived from the four foregroundcleaned CMB maps are in good agreement with those obtained by CamSpec and Plik using crossspectra; departures are well within 1σ for most parameters.
Fig. 12
Residuals of all mapbased bestfit models relative to CamSpec bestfit model (assuming a prior on τ) for ℓ_{max} = 2000. 
Inspecting the differences between the bestfit models derived from the four foregroundcleaned maps and from CamSpec plotted in Fig. 12, we find that the relative residuals are within 40 μK^{2} for all multipole ranges, and smaller than 20 μK^{2} at high ℓ. This can be compared to the corresponding residuals for the FFP6 simulation shown in Appendix E.
The likelihood used for this analysis does not take into account some systematic effects that will affect our foregroundcleaned CMB maps, such as relative calibration uncertainties between the frequency channel maps used to construct them, or their beam uncertainties. These effects are accounted for in the likelihoods in Planck Collaboration XV (2014). We have also adopted a very simple twoparameter model for the residual extragalactic foregrounds. Despite these limitations, the four CMB maps yield cosmological parameters in agreement with the crossspectrum based likelihoods for a basic sixparameter ΛCDM model. Thus we can be confident that the CMB maps are consistent with the power spectrum analysis.
Before concluding this section, we show in Fig. 13 the results from a twoparameter fit of an amplitudetilt model to each of the four foregroundcleaned maps, downgraded to 6° and repixelized at an N_{side} = 32 grid. Clearly, the maps are virtually identical on large angular scales measured relative to cosmic variance, with any differences being smaller than 0.1σ in terms of cosmological parameters. However, it is worth noting that the bestfit model, (q,n) = (1,0), is in some tension with the lowℓ spectrum, at about 1.7σ in this plot. The same tension between large and small angular scales is observed in Planck Collaboration XV (2014) and Planck Collaboration XVI (2014) with higher statistical significance using the full Planck likelihood. Irrespective of physical interpretation, the calculations presented here demonstrate that these lowℓ features are robust with respect to component separation techniques.
Fig. 13
Lowℓ power spectrum amplitude and tilt constraints measured relative to the bestfit PlanckΛCDM model derived from foregroundcleaned CMB maps smoothed to 6° FWHM. The cross shows the bestfit model (q,n) = (1,0). 
7. Higherorder statistics
The foregroundcleaned CMB maps presented in this paper are used as inputs for most Planck analyses of higherorder statistics, including nonGaussianity studies (Planck Collaboration XXIV 2014), studies of statistical isotropy (Planck Collaboration XXIII 2014), gravitational lensing by largescale structure (Planck Collaboration XVII 2014), and of the ISW effect (Planck Collaboration XXIV 2014). In this section we provide a summary of the nonGaussianity and gravitational lensing results.
7.1. NonGaussianity
Primordial nonGaussianity is typically constrained in terms of the amplitude, , of the quadratic corrections to the gravitational potential, as well as by means of the threepoint correlation function based on different triangle configurations. The results from these calculations for the foregroundcleaned CMB maps are presented in Planck Collaboration XXIV (2014). After subtraction of the lensingISW correlation contribution, the final result is , as estimated from the SMICA map using the KSW bispectrum estimator (Komatsu et al. 2005), consistent within 1σ with results from other methods and foregroundcleaned maps.
Uncertainties are evaluated by means of the FFP6 simulations, and potential biases are studied using both Gaussian and nonGaussian CMB realizations. In particular, when a detectable level of primordial nonGaussianity () is injected into the FFP6 simulations, each foregroundcleaned map yielded a positive detection within 2σ of the expected value, recovering values of , 19.0 ± 7.5, 11.1 ± 7.6 and 19.7 ± 7.4 for CommanderRuler, NILC, SEVEM, SMICA, respectively. We see that NILC and SMICA demonstrate the best recovery of the injected nonGaussianity, and we favoured the latter for nonGaussian studies for its faster performance over NILC. The foregroundcleaned CMB maps presented in this paper do not provide significant evidence of a nonzero value of , and realistic simulations show that the component separation methods do not suppress real nonGaussian signatures within expected uncertainties. The implications of these results in terms of early Universe physics are discussed in the relevant papers (Planck Collaboration XXIV 2014; Planck Collaboration XXII 2014).
7.2. Gravitational lensing by largescale structure
Gravitational lensing by the intervening matter imprints a nonGaussian signature in the CMB, which allows the reconstruction of the gravitational potential integrated along the line of sight to the last scattering surface. In Planck Collaboration XVII (2014), this effect has been detected at a high significance level (greater than 25σ) using the Planck temperature maps. Specifically, the lensing induced correlations between the total intensity and its gradients have been used to reconstruct a nearly full sky map of the lensing potential φ, which has been used for further studies on Planck data, including the detection of a nonzero correlation with the ISW (Planck Collaboration XXIV 2014; Planck Collaboration XIX 2014) and other tracers of largescale structure (notably, significant correlation with the CIB is reported in Planck Collaboration XVIII 2014), as well as the estimate of the power spectrum of the lensing potential and the associated likelihood. The latter was constructed using a simple minimum variance combination of the 143 and 217 GHz maps on about 70% of the sky, as well as subtracting dust contamination using the 857 GHzPlanck channel as a template (Planck Collaboration XVII 2014). These lensing results have improved the cosmological constraints from Planck (Planck Collaboration XVI 2014).
The foregroundcleaned CMB maps described in Sect. 5 were used to perform a lensing extraction on a larger sky fraction, reaching about 87% of the sky. We found the lensing power spectrum to be in good agreement with the one obtained using the minimum variance combination, i.e., the signal agrees within 1σ in the majority of the angular domain bins, and is characterized by an equivalent uncertainty. The foregroundcleaned maps were further exploited on the baseline 70% sky fraction for assessing the robustness of the main reconstruction against the foreground contamination (Planck Collaboration XVII 2014).
We show that the component separation algorithms presented in this paper do not bias the lensing reconstruction in the case of the large sky fraction considered here. We consider FFP6 simulations including noise and lensed CMB signal, propagated through each of the component separation algorithms described in Sect. 3. We perform a lensing potential reconstruction in the pixel domain based on the CMB maps processed by the four component separation methods using the metis algorithm described in Planck Collaboration XVII (2014). This method uses the quadratic estimator presented in Okamoto & Hu (2003), which corrects for the meanfield bias caused by extra sources of statistical anisotropy in addition to the CMB.
For each method, we combine the masks of CO regions, nearby galaxies and compact objects as defined in Planck Collaboration XVII (2014), with the CG90 mask described in Sect. 4. This procedure results in masks with sky fractions f_{sky} = 0.836,0.851,0.850,0.846 for CommanderRuler, NILC SEVEM, and SMICA, respectively.
We estimate the lensing potential power spectrum, , following the methodology described in Planck Collaboration XVII (2014). It consists of a pseudoC_{ℓ} estimate based on a highlyapodized version of the lensing potential reconstruction, which has an effective available sky fraction f_{sky,2} = 0.648,0.690,0.686,0.683 for CommanderRuler, NILC, SEVEM and SMICA, respectively. The bandpower reconstructions in 17 bins in the range 2 ≤ ℓ ≤ 1025 are plotted in Fig. 14, as well as the residuals relative to the theoretical lens power spectrum. All algorithms achieved an unbiased estimation of the underlying lensing power spectrum, with χ^{2} = 10.58,17.34,18.54,15.30, for CommanderRuler, NILC, SEVEM, and SMICA respectively, with 17 degrees of freedom. The associated probabilitytoexceed (PTE) values are 83%, 36%, 29%, 50%.
The power spectrum estimates are in remarkable agreement with each other. However, the CommanderRuler solution has significantly larger uncertainties, as expected from its lower signaltonoise ratio to lensing due to its larger beam. These results on simulated foregroundcleaned CMB maps demonstrate that the component separation algorithms do not alter the lensing signal, and this provides a strategy for achieving a robust lensing reconstruction on the largest possible sky coverage. The foregroundcleaned maps have been used in Planck Collaboration XVII (2014) to obtain lensing potential estimates on 87% of the sky.
8. Foreground components
In this section we consider the diffuse Galactic components, and present fullsky maps of thermal dust and CO emission, as well as a single lowfrequency component map representing the sum of synchrotron, AME, and freefree emission. Our allsky CO map is a “type 3” product as presented in Planck Collaboration XIII (2014). To assess the accuracy of these maps, we once again take advantage of the FFP6 simulation. The CommanderRuler method used in the following is described in Appendix A and consists of a standard parametric Bayesian MCMC analysis at low angular resolution, followed by a generalized leastsquares solution for component amplitudes at high resolution.
8.1. Data selection and processing
Fig. 14
Lensing power spectrum estimates from FFP6 simulations using an apodized mask covering f_{sky,2} ≃ 0.70 of the sky. 
Estimated monopoles and dipoles in Galactic coordinates, all measured in thermodynamic μK.
We only use the seven lowest Planck frequencies, from 30 to 353 GHz. The two highest channels have significantly different systematic properties than the lower frequency bands, for instance concerning calibration, ZLE, and noise correlations, and they are more relevant to thermal dust and CIB studies than to the present CMB analysis. Studies of specific foregrounds (CO, thermal dust, CIB etc.) using all Planck frequencies as well as ancillary data are discussed in companion papers Planck Collaboration VI (2014), Planck Collaboration XXI (2014), Planck Collaboration XI (2014), and additional future publications will consider extensions to AME, synchrotron and freefree emission.
In order to obtain unbiased estimates of the spectral parameters across all frequency bands, each map is downgraded from its native resolution to a common angular resolution of 40′ and repixelized at N_{side} = 256, a limit imposed by the LFI 30 GHz channel. Once the spectral indices have been determined, we reestimate the component amplitudes at native Planck resolution (see Appendix A).
Although the smoothing operation introduces noise correlations between pixels, we model the noise of the smoothed maps as uncorrelated white noise with an effective standard deviation, σ(p), for each pixel p. This approximation does not bias the final solution, because the analysis is performed independently for each pixel. However, it is important to note that correlations between pixels are not taken into account in this analysis. The effective noise uncertainty, σ(p), is estimated using realistic noise simulations downgraded in the same way as the data. The measured instrumental bandpasses are taken into account by integrating the emission laws over the bandpass for each component at each Monte Carlo step in the analysis.
The monopole (zeropoint) of each frequency map is not constrained by Planck, but is rather determined by postprocessing, and associated with a nonnegligible uncertainty (see Table 5 of Planck Collaboration I 2014). In addition, each frequency map includes a significant monopole contribution from isotropic extragalactic sources and CIB fluctuations not traced by local Galactic structure, ranging from less than about 10–20 μK at 70 GHz to several hundreds of μK at 353 GHz. Finally, the effective dipole in each map is associated with significant uncertainty due to the large kinematic CMB dipole. In order to prevent these effects from introducing modelling errors during component separation, they must be fit either prior to or jointly with the Galactic parameters. Unfortunately, when allowing free spectral parameters per pixel, there is a nearperfect degeneracy among the offsets, the foreground amplitudes and the spectral indices, and in order to break this degeneracy, it is necessary to reduce the number of spectral degrees of freedom.
We adopt the method described by Wehus et al. (in prep.) for this purpose, which has the additional advantage of making minimal assumptions about the foreground spectra. In short, this method uses linear regression between data from CMBsubtracted maps evaluated on pixels falling within each large N_{side} = 8 pixel to estimate the relative offsets, m_{1} and m_{2}, between any two maps at each position on the sky. Each regression provides a constraint of the form m_{1} = am_{2} + b, where a and b are the slope and offset, respectively, and where each value of m_{i} consists of the sum of both a monopole and a dipole term evaluated at that position. The individual monopoles and dipoles can then be reconstructed by measuring a and b in different regions of the sky, exploiting spatial variations in spectral indices, and solving jointly for two monopoles and dipoles, including constraints from all positions. To minimize degeneracies, a positivity prior is imposed on the fit, such that statistically significant negative pixels are heavily penalized. For 44 and 70 GHz, we retain the dipole values determined during the mapmaking process, and do not attempt to fit them.
The resulting complete set of monopole and dipole values is listed in Table 2. As a crosscheck, we performed a dedicated Commander run in which we fitted for the dipole at 353 GHz, together with the foreground amplitudes and spectral indices, and only found sub μK differences. This channel is by far the most problematic in our data set in terms of offset determination, because of the very bright dust emission at this frequency. As a result, there is a large relative uncertainty between the zerolevel of the dust amplitude map and the 353 GHz channel offset not accounted for in the following analyses. However, the sum of the two terms is well determined, and a potential error in either therefore does not compromise the quality of the other signal components (e.g., CMB and lowfrequency components). A similar comment applies between the offset at 30 GHz and the zerolevel of the lowfrequency component, although at a significantly lower level. These degeneracies can only be broken by including additional observations at lower and higher frequencies, respectively, and this will be addressed in a future publication.
8.2. Component models and priors
Our model for the lowresolution CMB analysis includes four independent physical components: CMB; “lowfrequency” emission; CO emission; and thermal dust emission. It can be written schematically in the form (1)where A_{i}(p) denotes the signal amplitude for component i at pixel p, ν_{0,i} is the reference frequency for each component, and ν refers to frequency. (Note that for readability, integration over bandpass, as well as unit conversions between antenna, flux density and thermodynamic units, is suppressed in this expression.) Thus, each component is modelled with a simple frequency spectrum parameterized in terms of an amplitude and a small set of free spectral parameters (a powerlaw index for the lowfrequency component, and an emissivity index and temperature for the thermal dust component); no spatial priors are imposed. One goal of the present analysis is to understand how well this simple model captures the sky signal in terms of effective components over the considered frequency range, and we exploit the FFP6 simulation (see Sect. 4) for this purpose.
In order to take into account the effect of bandpass integration, each term in the above model is evaluated as an integral over the bandpass as described in Sect. 3 of Planck Collaboration IX (2014), and converted internally to thermodynamic units. Accordingly, the reference frequencies in Eq. (1) are computed as effective integrals over the bandpass, such that the amplitude map, A_{i}(p), corresponds to the foreground map observed by the reference detector, i.e., after taking into account the bandpass. In order to minimize degeneracies between the different signal components, the reference band for a given component is set to the frequency at which its relative signaltonoise ratio is maximized.
The foreground model defined in Eq. (1) is motivated by prior knowledge about the foreground composition over the CMB frequencies as outlined in Sect. 2, as is our choice of priors. In addition to the Jeffreys prior^{5} (Eriksen et al. 2008), we adopt Gaussian priors on all spectral parameters with centre values and widths attempting to strike a balance between prior knowledge and allowing the data to find the optimal solution. Where needed, we have also run dedicated analyses, either including particular high signaltonoise ratio subsets of the data or using a lower resolution parameterization to increase the effective signaltonoise in order to inform our prior choices. We now consider each foreground component in turn, and note in passing that the CMB component, by virtue of being a blackbody signal, is given by a constant in thermodynamic temperature units.
We approximate the lowfrequency component by a straight power law in antenna temperature with a free spectral index per pixel, and adopt a prior of β = −3 ± 0.3 (this is the index in terms of antenna temperature). This choice is determined by noting that the prior is in practice only relevant at high Galactic latitudes where the signaltonoise ratio is low and the dominant foreground component is expected to be synchrotron emission; in the signaldominated and lowlatitude AME and freefree regions, the data are sufficiently strong to render the prior irrelevant. For validation purposes, we have also considered minor variations around this prior, such as β = −2.9 ± 0.3 and β = −3.05 ± 0.2, finding only small differences in the final solutions. The reference band for the lowfrequency component is set to 30 GHz, where the lowfrequency foreground signal peaks. The final lowfrequency amplitude map is provided in units of thermodynamic microkelvin.
Fig. 15
Posterior mean foreground amplitude maps derived from the lowresolution analysis. From top to bottom are shown the lowfrequency, CO and thermal dust emission maps. 
Fig. 16
Posterior mean spectral parameter maps derived from the lowresolution analysis. The top panel shows the power law index of the lowfrequency component, and the bottom panel shows the emissivity index of the onecomponent thermal dust model. Note that the systematic error due to monopole and dipole uncertainties is significant for the dust emissivity in regions with a low thermal dust amplitude. 
Fig. 17
χ^{2} per pixel for the joint CMB and foreground analysis. The expected value for an acceptable fit is 7, corresponding to the number of frequency bands used in this analysis. The pixels with high values can be classified into two types, due to either modelling errors (i.e., high residuals in the Galactic plane) or to unmodelled correlated noise (i.e., stripes crossing through low dust emission regions). 
The CO emission is modelled in terms of a single line ratio for each frequency. Specifically, the CO amplitude is normalized to the 100 GHz band, and defined in units of K km s^{1} (Planck Collaboration XIII 2014). The amplitude at other frequencies is determined by a single multiplicative factor relative to this, with a numerical value of 0.595 at 217 GHz and 0.297 at 353 GHz; all other frequencies are set to zero. These values are obtained from a dedicated CO analysis that includes only high signaltonoise ratio CO regions covering a total of 0.5% of the sky. The derived values are in good agreement with those presented by Planck Collaboration XIII (2014).
Thermal dust emission is modelled by a onecomponent modified blackbody emission law with a free emissivity spectral index, β_{d}, and dust temperature, T_{d}, per pixel. However, since we only include frequencies below 353 GHz, the dust temperature is largely unconstrained in our fits, and we therefore impose a tight prior around the commonly accepted mean value of T_{d} = 18 ± 0.05 K. The only reason we do not fix it completely to 18 K is to allow for modelling errors near the Galactic centre. The dust emissivity prior is set to β_{d} = 1.6 ± 0.3, where the mean is determined by a dedicated run fitting for a single bestfit value for the highlatitude sky, where the prior is relevant. The reference band for the thermal dust component is 353 GHz, and the final map is provided in units of megajansky per steradian.
8.3. Results and validation
The output of the Bayesian component separation algorithm is a set of samples drawn from the joint posterior distribution of the model parameters, as opposed to a single welldefined value for each. For convenience, we summarize this distribution in terms of posterior mean and standard deviation maps, computed over the sample set, after rejecting a short burnin phase. The goodnessoffit is monitored in terms of the χ^{2} per pixel. Although convenient, it is, however, important to note that this description does not provide a comprehensive statistical representation of the full posterior distribution, which is intrinsically nonGaussian. One should be careful about making inferences in the low signaltonoise regime based on this simplified description.
The lowresolution Commander posterior mean amplitude maps are shown in Fig. 15 for the lowfrequency, CO, and thermal dust components, and the spectral index maps in Fig. 16. The associated χ^{2} map is plotted in Fig. 17. Note that because we are sampling from the posterior instead of searching for the maximumlikelihood point, the expected number of degrees of freedom is equal to N_{band} = 7 in this plot, not N_{band} − N_{par}.
Several features can be seen here, foremost of which is that the Galactic plane is strikingly obvious, with χ^{2} values exceeding 10^{4} for seven degreesoffreedom in a few pixels. This is not surprising, given the very simplified model at low frequencies (i.e., a single power law accounting for AME, synchrotron, and freefree emission), as well as the assumption of a nearly constant dust temperature of 18 K. Second, there is an extended region of moderately high χ^{2} roughly aligned with great circles going through the ecliptic poles, indicating the presence of correlated noise in the scanning rings not accounted for in our white noise model.
Based on these, and other considerations, it is clear that parts of the sky must be masked before proceeding to CMB power spectrum and likelihood analyses. This masking process is discussed at greater length in Planck Collaboration XV (2014), and results in different masks for specific applications. The goal of our present discussion is to evaluate the adequacy of the mask adopted for lowℓ likelihood analysis (L87), which is based on the fits presented here. This mask removes 13% of the sky, and is derived from a combination of χ^{2} and component amplitude thresholding.
Figure 18 compares the highresolution Ruler solution to the lowresolution Commander solution for CMB, CO and thermal dust on a particularly strong CO complex near the Fan region, centred on Galactic coordinates (l,b) = (110°,15°).
For validation purposes, we analyse the simulations described in Sect. 4 in the same way as the real data, including monopole and dipole determination, CO line ratio estimation and spectral index estimation. Individual component maps at each observed frequency are available from the simulation process, and used for direct comparison with the reconstructed products.
In Fig. 19 we show the differences between the recovered and input component maps at their respective reference frequencies. The boundary of the 13% Commander mask is traced by the white contours, and a bestfit monopole and dipole have been subtracted from each difference map. All difference maps are shown in units of thermodynamic μK. The top panel of Fig. 20 gives the error histograms outside the masked region for each component, normalized to the respective estimated standard deviation; if the recovered solution has both correct mean and standard deviation, these histograms should match a Gaussian distribution with zero mean and unit variance, indicated by the dashed black line. Conversely, a significant bias would be visible as a horizontal shift in this plot, while underestimation of the errors would result in too wide a distribution and vice versa. The bottom panel shows the fractional error (i.e., the error divided by the true input value) for all pixels with signal above 5σ; the fractional error is not a useful quantity for noisier signals.
Fig. 18
Comparison of the highresolution Ruler (top) and lowresolution Commander (bottom) amplitude maps for a particularly strong CO complex near the Fan region; the maps are centred on Galactic coordinates (l,b) = (110°,15°), and the grid spacing is 5°. Columns show, left to right: the CMB amplitude; the CO amplitude at 100 GHzand the thermal dust amplitude at 353 GHz. 
The difference maps in Fig. 19 display significant errors in the Galactic plane. For the lowfrequency component, the residuals are dominated by freefree emission, while for thermal dust the dominant contaminant is CO emission. However, outside the mask the residuals are small, and, at least for the lowfrequency and CO components, the spatial characteristics appear similar to instrumental noise. This is more clear in the histograms shown in the top panel of Fig. 20; the mean and standard deviations are δ_{lf} = 0.01 ± 1.12, δ_{CO} = 0.00 ± 0.87, and δ_{td} = 0.00 ± 2.01, respectively, for the lowfrequency, CO and thermal dust components. There is no evidence of bias outside the mask in any component, and the error estimates are accurate to 12 and 13% for the lowfrequency and CO components. Note, though, that the estimated error for the CO component is actually larger than the true uncertainty, suggesting that the white noise approximation for the 100 GHz channel overestimates the true noise. This can occur if the correlated instrumental noise is important in regions where there is no significant CO emission. Locally rescaling the white noise to account for spatially varying correlated noise would correct this effect.
For the thermal dust component, on the other hand, the error is underestimated by a factor of 2. The explanation for this is most easily seen from the lower panel of Fig. 19. This map is dominated by isotropic CIB fluctuations, rather than instrumental noise. Because these fluctuations have a slightly different spectrum than the dominant Galactic dust emission, and the model does not account for a separate CIB component, the error on the Galactic component is underestimated. When using the Galactic map presented here for detailed analysis near the noise limit, taking into account these residual fluctuations is essential, and the effective noise per pixel should be increased by a factor of 2.
As clearly seen in Fig. 19, the residuals inside the mask are highly significant in a strict statistical sense. However, as seen in the bottom panel of Fig. 20, they are relatively small in terms of fractional errors. Specifically, the three histograms have means and standard deviations of f_{lf} = 0.00 ± 0.10, f_{CO} = −0.03 ± 0.10, and f_{td} = 0.00 ± 0.06, respectively, for the lowfrequency, CO and thermal dust components. The largest bias is observed for the CO component, for which the absolute amplitude is biased by 3%. The bias in the lowfrequency and thermal dust components is negligible, and the fractional uncertainties are 10 and 6%, respectively. This confirms that approximating the sum of the three lowfrequency components by a single powerlaw over the Planck frequency bands is reasonable; if modelling errors dominated, one would expect to see a significant bias in the resulting amplitude.
In order to validate the spectral parameters, we show in Fig. 21 histograms of the normalized residuals for each foreground component evaluated at its two leading subdominant frequencies (i.e., at 44 and 70 GHz for the lowfrequency component; at 217 and 353 GHz for the CO component; and at 143 and 217 GHz for the thermal dust component). The means and standard deviations of these distributions are: δ_{lf}(44 GHz) = −0.41 ± 1.98 and δ_{lf}(70 GHz) = −0.34 ± 2.04 for the lowfrequency component; δ_{CO}(217 GHz) = 0.10 ± 0.84 and δ_{CO}(353 GHz) = 0.51 ± 1.00 for the CO component; and δ_{td}(143 GHz) = −0.02 ± 1.53 and δ_{td}(217 GHz) = −0.13 ± 1.87 for the thermal dust component. As expected, the effect of modelling errors is more significant at the subdominant frequencies than at the pivot frequencies, when measured in terms of statistical uncertainties, since the foreground signal is weaker and the confusion with the other components relatively larger. Nevertheless, we see that the absolute bias is at most 0.5σ for the CO component at 353 GHz, while the thermal dust bias is negligible even at 143 GHz. The estimated uncertainties are generally underestimated by up to a factor of two due to these modelling errors.
Fig. 19
Amplitude residual maps, A_{out} − A_{in}, computed blindly from the FFP6 simulation. The panels show (from top to bottom) the lowfrequency residual at 30 GHz, the CO residual at 100 GHz and the thermal dust residual at 353 GHz. All units are thermodynamic μK. The white lines indicate the boundary of the Commander likelihood analysis mask, removing 13% of the sky. 
Fig. 20
Error validation for component amplitudes, evaluated from the FFP6 simulation. The upper panel shows histograms of the normalized errors δ = (A_{out} − A_{in}) /σ_{out} for the three foreground components and including all pixels outside the Commander likelihood analysis mask. The lower panel shows histograms of the fractional error f ≡ (A_{out} − A_{in}) /A_{in} for pixels with a foreground detection level above 5σ. No evidence of significant bias is observed for any component, and the uncertainty estimates for the lowfrequency and CO components are accurate to about 12%; the thermal dust uncertainty is underestimated by a factor of 2 due to the presence of unmodelled fluctuations. 
Fig. 21
Validation of spectral parameters for lowfrequency foregrounds, thermal dust, and CO emission, evaluated from the FFP6 simulation. Each histogram shows the error distribution at the two leading subdominant frequencies in the form of the normalized errors δ = (A_{out}(ν) − A_{in}(ν)) /σ_{out}(ν) for all pixels outside the Commander likelihood analysis mask, where A_{out}(ν) is the predicted foreground amplitude at frequency ν given the estimated amplitude and spectral parameters, and σ_{out}(ν) is the corresponding standard deviation computed over the sample set. 
Finally, the efficiency of the adopted foreground model for CMB analysis is quantified in Appendix E in terms of power spectrum residuals and cosmological parameter estimation.
To summarize, we find that the simplified model, defined by Eq. (1), provides a good fit to the realistic FFP6 simulation for most of the sky. Absolute residuals are small, and the amplitude uncertainty estimates are accurate to around 12%, except for the thermal dust component for which unmodelled CIB fluctuations are important. Further, we find that the real Planck data behave both qualitatively and quantitatively very similarly to the FFP6 simulation, suggesting that this approach also performs well on the real sky.
8.4. Interpretation and comparison with other results
The maps shown in Figs. 15 and 16 provide a succinct summary of the average foreground properties over the Planck frequency range. We now consider their physical interpretation and compare them to products from alternative methods.
First, the top panel of Fig. 22 shows a difference map between the thermal dust map^{6} derived in the present paper using Planck frequencies between 30 and 353 GHz, and the one determined from only the three highest Planck frequencies (353 to 857 GHz) and the 100 μm IRIS map by Planck Collaboration XI (2014). The largescale features in this map are dominated by the different CO, zodiacal light and offset modelling approaches adopted by the two pipelines. Specifically, while the Commander solution jointly fits for CO and thermal dust, the highfrequency solution assumes CO to be negligible at these frequencies; while Commander analyzes raw Planck maps with ZLE still present, the highfrequency analysis considers maps that have been explicitly corrected for zodiacal light using the Planck Collaboration VI 2014 model (which itself has a poorly defined zerolevel); and while the Commander solution fits internally for spurious monopoles and dipoles by regression as described in Sect. 8.1, the highfrequency solution determines these through correlations with external HI observations. The middle panel of Fig. 22 shows the difference between the two solutions after subtracting residual monopoles, dipoles and zodiacal light components. The differences are now below 0.01 MJy sr^{1} everywhere at high latitudes. This is to be compared with the corresponding thermal dust amplitude map in the bottom panel of Fig. 15, which ranges between 0 and 2.5 MJy sr^{1}. Inside the Galactic plane, the differences are dominated by residuals due to different CO modelling, seen as solid blue colours in Fig. 22; however, even in this region the differences are smaller than 5% of the amplitude. Thus, the agreement between the two solutions is excellent, despite their very different choice of data and processing methodology. The bottom panel shows the T–T plot between the two maps, excluding any pixel for which the Commander CO amplitude at 100 GHz is larger than 1 K km s^{1}.
Fig. 22
Top: difference map between the thermal dust amplitude at 353 GHz presented in this paper using the Planck 30 to 353 GHz frequencies, and that derived by Planck Collaboration XI (2014) using only the Planck 353, 545 and 857 GHz channels and the 100 μm IRIS map. Both maps are smoothed to a common resolution of 40′. Middle: the same difference map, but accounting for relative monopole, dipole and zodiacal emission treatment. Bottom: T–T plot between the same two maps after applying relative corrections. Any pixels with a CO amplitude at 100 GHz larger than 1 K km s^{1} are removed from this plot. 
From Fig. 16, we see that the dust emissivity ranges between 1.3 and 1.7 for most of the sky, in good agreement with the mean value of β_{d} = 1.56 ± 0.02 derived by Planck Collaboration Int. XXII (2014) by averaging over 10° disks and including additional ancillary data. Considering only the pixels with a posterior distribution width that is a third of the prior width (i.e., σ(β_{d}) < 0.1), we find a mean value of 1.49. The two exceptions are a large region of shallow indices northeast of the Galactic centre, and steep indices near the Galactic plane. The former region corresponds to a part of the sky with low dust emission, where we expect the spectral index to be sensitive to both monopole and dipole residuals, as well as instrumental systematics, such as correlated 1 /f noise. The latter appears to be particularly pertinent here because the shallow index region at least partially traces the Planck scanning strategy; as a result, the systematic error on the spectral index in this region is considerable. The main systematic uncertainty connected to the region of steep indices around the Galactic plane is confusion with CO emission.
The CO map shown in Figs. 15 and 18 is discussed in greater detail in Planck Collaboration XIII (2014). A distinct advantage of this particular solution over available alternatives is its high signaltonoise ratio per pixel, which is achieved by reducing all information into a single value per pixel. Consequently, this map serves as a unique tool for followup CO observations. However, the assumption of a constant line ratio over the full sky may lead to a significant systematic uncertainty on CO amplitude per pixel.
Finally, the spectral index map for the lowfrequency component shown in Fig. 16 can be used to determine the dominant lowfrequency component (synchrotron, freefree or AME) as a function of position on the sky. To illustrate this connection, we once again take advantage of the FFP6 simulation for which we know the amplitude of each lowfrequency component per pixel. In the top panel of Fig. 23, we use this information to make a “dominant component map”; dark blue indicates that synchrotron emission is strongest at a given pixel, light blue that freefree is strongest, and orange that spinning dust (AME) dominates. In the bottom panel, we show our derived powerlaw index map from the same simulation. As expected, the correspondence between the powerlaw index and the dominant component is very strong, implying that the spectral index map can be used to trace the individual components. In particular, we see that an index below about − 3.3 reflects the presence of a component consistent with a spinning dust model peaking below 20 GHz over the Planck frequency range, while an index higher than around − 2.3 signals the importance of freefree emission. Intermediate values typically indicate synchrotron emission, although it should be noted that the signaltonoise ratio at very high latitudes is low and the results are therefore priordriven in these regions.
Returning to the spectral index map shown in Fig. 16, we see a good correspondence between the real data and the simulation. Features present in the simulation also appear in the data. For instance, we see that the spectral index in the socalled Fan Region (i.e., near Galactic coordinates (l,b) = (90°,20°)) is low in both cases, and this alone provides strong evidence for the presence of AME. Further, the AME spectral index is consistent with the spinning dust interpretation. This powerlaw index map may be used to identify particular AME regions for followup observations. Finally, we note that regions known for strong freefree emission, such as the Gum Nebula or Zeta Ophiuchi, have spectral indices close to − 2.1 or − 2.2, as expected.
9. Conclusions
Fig. 23
Top: dominant foreground component per pixel at 30 GHz in the FFP6 simulation. Dark blue indicates that synchrotron emission is the strongest component at 30 GHz, light blue indicates that freefree dominates, and orange indicates that spinning dust (AME) is the strongest component. Bottom: the recovered lowfrequency powerlaw index derived from the same simulation. 
We have tested four component separation algorithms on the Planck frequency maps to produce clean maps of the CMB anisotropies over a large area of sky. These CMB maps are used for studies of statistics and isotropy (Planck Collaboration XXIII 2014), primordial nonGaussianity (Planck Collaboration XXIV 2014), gravitational lensing (Planck Collaboration XVII 2014), the ISW effect (Planck Collaboration XIX 2014), cosmic geometry and topology (Planck Collaboration XXVI 2014), searches for cosmic defects from primordial phase transitions (Planck Collaboration XXV 2014), as well as an integral part of the lowℓPlanck likelihood (Planck Collaboration XV 2014). Two of the methods, one using internal foreground templates (SEVEM) and the other an ILC in needlet space (NILC), are nonparametric, extracting the CMB map by minimizing the variance of the total contamination. The other two methods fit models of the foregrounds to clean the CMB of their emission. One fits a parametric model in real space (CR) and one fits a nonparametric in the harmonic domain (SMICA).
All four methods have been demonstrated to work well both on real and simulated data, and to yield consistent results. Nevertheless, there are differences between the methods, making them more or less suitable for specific applications. For instance, CommanderRuler allows a joint parametric foreground estimation and CMB power spectrum estimation, with full propagation of foreground uncertainties to cosmological parameters, but is limited to a lower angular resolution than the other codes. This method has therefore been selected for the lowℓPlanck likelihood (Planck Collaboration XV 2014) and to produce astrophysical component maps (Sect. 8), while it is suboptimal for applications requiring full angular resolution, e.g., gravitational lensing reconstruction or estimation of primordial nonGaussianity. For these purposes, we use the three higherresolution maps. We take SMICA to be the leading method, based on its superior performance on the FFP6 simulation, where it has be shown to have the lowest residual foreground contamination at large scales and to preserve primordial nonGaussianity. When subjecting foregroundcleaned Planck maps to scientific analysis, we use the other two or three maps, as appropriate, to assess the uncertainties inherent in the choice of methods and the assumptions they make. Indeed, this is the main purpose for presenting four different CMB solutions to the general community.
The CMB anisotropies are robustly recovered over a large fraction (73%) of the sky and down to small angular scales, reaching to multipoles ℓ ≈ 2000. We characterize the CMB maps with angular power spectra and cosmological parameter constraints. Parameter constraints from these maps are consistent with those from the Planck likelihood function based on crossspectra and large sky cuts (Planck Collaboration XV 2014). This agreement supports the robustness of both our component separation methodology and cosmological parameter constraints.
The realspace parametric fits of CommanderRuler enable us to characterize the diffuse Galactic foregrounds. We parameterize them with a lowfrequency powerlaw component, representing the sum of synchrotron, freefree, and AME emission, a highfrequency modified blackbody spectrum describing thermal dust emission, and a molecular CO component. Using only the Planck data from 30 to 353 GHz, we fit for the amplitude and spectral parameters of the three foregrounds and the CMB simultaneously at each pixel of a 40arcmin resolution map. The spectral parameters are the lowfrequency component powerlaw exponent and the modified blackbody emissivity powerlaw exponent; the CO line ratios are spatially fixed. These parameters give us the source mixing matrix, which we then use in a direct inversion to deduce the component amplitudes at higher resolution. Through Gibbs sampling, we obtain realizations drawn from the full posterior distribution of possible foreground and CMB solutions, giving us a powerful ability to statistically characterize our results.
Our indepth analysis of the recovered CMB anisotropies is unprecedented for component separation studies, concerning both the accuracy of cosmological parameter constraints, and studies of early Universe physics and structure formation through gravitational lensing. On the other hand, the complex nature of the foreground emission over such a large frequency range limits us to the use of relatively simple methods when analysing Planck data alone. An extensive study in combination with other probes of Galactic foregrounds will be presented in forthcoming papers. In particular, the separation of individual components at low frequencies requires the use of ancillary data, for example, from the WMAP and radio surveys.
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 (in particular the lead countries France and Italy), with contributions from NASA (USA) and telescope reflectors provided by a collaboration between ESA and a scientific consortium led and funded by Denmark.
Note that SMICA, being defined in harmonic space, employs a smooth filling process inside a small Galactic mask to prevent foreground residuals from leaking from low to high Galactic latitudes, and therefore appears visually different from the other three solutions in this respect; see Appendix D.
Note that while the Commander component maps are normalized with respect to effective Planck frequencies (i.e., including bandpass integration; see Sect. 8.2), the highfrequency dust map is normalized to nominal Planck frequencies. In this comparison the former has been multiplied by 0.987 to account for this difference.
Acknowledgments
The development of Planck has been supported by: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. The authors acknowledge the support provided by the Advanced Computing and eScience team at IFCA. This work made use of the COSMOS supercomputer, part of the STFC DiRAC HPC Facility. Some of the results in this paper have been derived using the HEALPix package.
References
 AliHaïmoud, Y., Hirata, C. M., & Dickinson, C. 2009, MNRAS, 395, 1055 [NASA ADS] [CrossRef] [Google Scholar]
 Banday, A. J., Dickinson, C., Davies, R. D., Davis, R. J., & Górski, K. M. 2003, MNRAS, 345, 897 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97 [NASA ADS] [CrossRef] [Google Scholar]
 BenoitLévy, A., Déchelette, T., Benabed, K., et al. 2013, A&A, 555, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Betoule, M., Pierpaoli, E., Delabrouille, J., Le Jeune, M., & Cardoso, J.F. 2009, A&A, 503, 691 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonaldi, A., Ricciardi, S., Leach, S., et al. 2007, MNRAS, 382, 1791 [NASA ADS] [Google Scholar]
 Broadbent, A., Osborne, J. L., & Haslam, C. G. T. 1989, MNRAS, 237, 381 [NASA ADS] [CrossRef] [Google Scholar]
 Cardoso, J.F., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, IEEE J. Select. Topics Signal Process., 2, 735 [Google Scholar]
 Casaponsa, B., Barreiro, R. B., Curto, A., MartínezGonzález, E., & Vielva, P. 2011, MNRAS, 411, 2019 [NASA ADS] [CrossRef] [Google Scholar]
 Davies, R. D., Watson, R. A., & Gutierrez, C. M. 1996, MNRAS, 278, 925 [NASA ADS] [CrossRef] [Google Scholar]
 Davies, R. D., Dickinson, C., Banday, A. J., et al. 2006, MNRAS, 370, 1125 [NASA ADS] [CrossRef] [Google Scholar]
 de OliveiraCosta, A., Tegmark, M., Davies, R. D., et al. 2004, ApJ, 606, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dickinson, C., Davies, R. D., & Davis, R. J. 2003, MNRAS, 341, 369 [Google Scholar]
 Dobler, G., & Finkbeiner, D. P. 2008, ApJ, 680, 1222 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T., & Lazarian, A. 1998, ApJ, 508, 157 [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]
 Eriksen, H. K., Dickinson, C., Lawrence, C. R., et al. 2006, ApJ, 641, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Huey, G., Saha, R., et al. 2007, ApJ, 656, 641 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Fendt, W., & Wandelt, B. 2008, in APS April Meeting Abstracts, J8004 [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. 2004, ApJ, 614, 186 [NASA ADS] [CrossRef] [Google Scholar]
 Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 1999, ApJ, 524, 867 [Google Scholar]
 Finkbeiner, D. P., Langston, G. I., & Minter, A. H. 2004, ApJ, 617, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Ghosh, T., Banday, A. J., Jaffe, T., et al. 2012, MNRAS, 422, 3617 [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., Stoffel, H., Salter, C. J., & Wilson, W. E. 1982, A&ASS, 47, 1 [NASA ADS] [Google Scholar]
 Hivon, E., Górski, K. M., Netterfield, C. B., et al. 2002, ApJ, 567, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Hoang, T., & Lazarian, A. 2012, Adv. Astron., 2012, 208159 [NASA ADS] [CrossRef] [Google Scholar]
 Jewell, J., Levin, S., & Anderson, C. H. 2004, ApJ, 609, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Kogut, A. 1996, in BAAS, 28, 1295 [Google Scholar]
 Komatsu, E., Spergel, D. N., & Wandelt, B. D. 2005, ApJ, 634, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Lagache, G. 2003, A&A, 405, 813 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leach, S. M., Cardoso, J., Baccigalupi, C., et al. 2008, A&A, 491, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leitch, E. M., Readhead, A. C. S., Pearson, T. J., & Myers, S. T. 1997, ApJ, 486, L23 [NASA ADS] [CrossRef] [Google Scholar]
 MartınezGonzalez, E., Diego, J. M., Vielva, P., & Silk, J. 2003, MNRAS, 345 [Google Scholar]
 MivilleDeschênes, M.A., Ysard, N., Lavabre, A., et al. 2008, A&A, 490, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Okamoto, T., & Hu, W. 2003, Phys. Rev., D, 67, 083002 [Google Scholar]
 Patanchon, G., Cardoso, J.F., Delabrouille, J., & Vielva, P. 2005, MNRAS, 364, 1185 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Górski, K. M., Bartlett, J., et al. 2012, ApJ, 755, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration XIII. 2011, A&A, 536, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2011, A&A, 536, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. IX. 2013, A&A, 554, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XII. 2013, A&A, 557, A53 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration 2013, The Explanatory Supplement to the Planck 2013 results, http://www.sciops.esa.int/wikiSI/planckpla/index.php?title=Main_Page (ESA) [Google Scholar]
 Planck Collaboration I. 2014, A&A, 571, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2014, A&A, 571, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2014, A&A, 571, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2014, A&A, 571, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2014, A&A, 571, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2014, A&A, 571, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2014, A&A, 571, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2014, A&A, 571, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [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 XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2014, A&A, 571, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2014, A&A, 571, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2014, A&A, 571, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXI. 2014, A&A, 571, A31 [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. XXII. 2014, A&A, submitted [Google Scholar]
 Platania, P., Burigana, C., Maino, D., et al. 2003, A&A, 410, 847 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reich, P., & Reich, W. 1988, A&AS, 74, 7 [NASA ADS] [Google Scholar]
 Rocha, G., Contaldi, C. R., Colombo, L. P. L., et al. 2010 [arXiv:1008.4948] [Google Scholar]
 Rocha, G., Contaldi, C. R., Bond, J. R., & Górski, K. M. 2011, MNRAS, 414, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Stompor, R., Leach, S., Stivoli, F., & Baccigalupi, C. 2009, MNRAS, 392, 216 [Google Scholar]
 Tegmark, M., de OliveiraCosta, A., & Hamilton, A. 2003, Phys. Rev. D, 68, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Tristram, M., Patanchon, G., MacíasPérez, J. F., et al. 2005, A&A, 436, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wandelt, B. D., Larson, D. L., & Lakshminarayanan, A. 2004, Phys. Rev. D, 70, 083511 [NASA ADS] [CrossRef] [Google Scholar]
 Ysard, N., & Verstraete, L. 2010, A&A, 509, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ysard, N., MivilleDeschênes, M. A., & Verstraete, L. 2010, A&A, 509, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Physical parametrization
The CommanderRuler (CR) approach implements Bayesian component separation in pixel space, fitting a parametric model to the data by sampling the posterior distribution for the model parameters. For computational reasons, the fit is performed in a twostep procedure: first, both foreground amplitudes and spectral parameters are found at lowresolution using Markov Chain Monte Carlo (MCMC)/Gibbs sampling algorithms (Jewell et al. 2004; Wandelt et al. 2004; Eriksen et al. 2004, 2007, 2008); second, the amplitudes are recalculated at high resolution by solving the generalized least squares system (GLSS) per pixel with the spectral parameters fixed to their values from the lowresolution run.
For the CMBoriented analysis presented in this paper, we only use the seven lowest Planck frequencies, i.e., from 30 to 353 GHz. We first downgrade each frequency map from its native angular resolution to a common resolution of 40 arcminutes and repixelize at HEALPixN_{side} = 256. Second, we set the monopoles and dipoles for each frequency band as described in Sect. 8.1, using a method that locally records spectral indices (Wehus et al., in prep.). We approximate the effective instrumental noise as white with a root mean square (rms) per pixel given by the Planck scanning pattern and an amplitude calibrated by smoothing simulations of the instrumental noise, including correlations, to the same resolution. For the highresolution analysis, the important preprocessing step is the upgrading of the effective lowresolution mixing matrices to full Planck resolution: this is done by repixelizing from N_{side} = 256 to 2048 in harmonic space, ensuring that potential pixelization effects from the lowresolution map do not introduce sharp boundaries in the highresolution map.
Our model for the data, a map d_{ν} at frequency ν, consists of a linear combination of N_{c} astrophysical components, plus instrumental noise, (A.1)where A^{i} denotes a sky map vector containing the foreground amplitude map for component i normalized at a reference frequency, and is a diagonal matrix describing the spectral emission law for component i as a function of frequency and which depends on a (small) set of spectral parameters, θ. The CMB signal is included in the sum and, as a special case, it may be represented either in harmonic or pixel space, depending on whether the main goal of the analysis is CMB power spectrum analysis or component separation. The former representation is used for the Commanderbased lowℓ likelihood presented in Planck Collaboration XV (2014), while the latter is used for the foreground fits presented in Sect. 8 of this paper.
Bayes theorem specifies the posterior distribution for the model parameters, (A.2)where ℒ(A^{i},θ) = P(A^{i},θ  d) is a Gaussian likelihood of observing data, d, given model parameters (A^{i},θ), and the prior P(A^{i},θ^{i}) depends on the application. In this paper, the prior on spectral indices is a product of a Jeffreys prior and physical priors, as detailed in Sect. 8.2; no priors are imposed on the foreground amplitudes.
In the lowresolution Commander analysis, we exploit a Gibbs sampler to map out the posterior distribution (Eriksen et al. 2008), adopting the following minimal twostep scheme: The first conditional distribution is a multivariate Gaussian, while the second distribution does not have an analytic form and must be mapped out numerically.
The highresolution Ruler analysis maximizes the foreground amplitude conditional in Eq. (A.3) numerically by solving the generalized least squares system (A.5)were N_{ν} is the noise covariance matrix (assumed to be diagonal) of the νth channel, F^{i} ≡ F^{i}(θ) and we have neglected the different angular resolutions of the channels. The posterior marginal average for the highresolution amplitude maps is then given by , a sum over the N_{sample} samples of the spectral parameters θ.
Once the channel weights, W_{ν}, have been computed, processing a large number of simulations requires negligible computational resources. This feature has been extensively used for computation of the effective beam of Ruler maps: FFP6 CMB simulations for the 30 to 353 GHz channels are combined according to W_{ν} and the effective beam transfer function is found as . Here, is the power spectrum of the input simulation before convolution with the instrumental beam, w_{ℓ} is the HEALPix pixel window function, and the average is taken over the set of simulations. Missing pixels are set to 0 when computing C_{ℓ} in the above expression. The low number (~500) of missing pixels in the data renders the impact of such a choice negligible at ℓ< 2000. A similar procedure is used for defining the effective beam of the nonCMB components.
The above algorithm produces a set of samples drawn from the posterior distribution, as opposed to a direct estimate of individual component amplitudes or spectral parameters. While this sample set provides a statistically complete representation of the posterior, it is nontrivial to visualize or to compare the distribution with external data. For convenience, we therefore summarize the distribution in terms of mean and standard deviation maps for each component. We emphasize, however, that the distribution is significantly nonGaussian, and when searching for features in the maps at low signaltonoise levels, one must take into account the exact distribution.
Finally, the CommanderRuler confidence mask (see Sect. 5) is primarily defined by the product of the CG80 mask and the point source mask described in Sect. 4. We additionally remove any pixels excluded by the 13% Commander likelihood mask described by Planck Collaboration XV (2014); however, this is almost entirely included within the CG80 mask, and this step therefore has very little impact on the final result. To complete the mask, we remove any pixels for which the highresolution Ruler CMB map, smoothed to 40 arcmin, differs by more than 3σ from the lowresolution Commander CMB map, which can happen due to spatial spectral variations on pixel scales.
Appendix B: Internal linear combination
NILC is a method to extract the CMB (or any component with known spectral behaviour) by applying the ILC technique to multichannel observations in needlet space, that is, with weights that are allowed to vary over the sky and over the full multipole range.
The ability to linearly combine input maps varying over the sky and over multipoles is called localization. In the needlet framework, harmonic localization is achieved using a set of bandpass filters defining a series of scales, and spatial localization is achieved at each scale by defining zones over the sky. The harmonic localization adopted here uses nine spectral bands covering multipoles up to ℓ = 3200 (see Fig. B.1). The spatial localization depends on the scale. At the coarsest scale, which includes the multipoles of lowest degree, we use a single zone (no localization), while at the finest scales (which include the highest multipoles) the sky is partitioned into 20 zones (again, see Fig. B.1).
Fig. B.1
Spectral localization for NILC using nine spectral window functions defining nine needlet scales (top panel). The scaledependent spatial localization partitions the sky into one zone (for scale 1), two zones (for scale 2), four zones (for scale 3), or twenty zones (for scales 5, 6, 7, 8 and 9). The twozone, 4zone and 20zone partitions are depicted in the lower panels. 
The NILC method amounts to computing an ILC in each zone of each scale, allowing the ILC weights to adapt naturally to the varying strength of the other components as a function of position and multipole. A complete description of the basic NILC method can be found in Delabrouille et al. (2009).
In the present work, however, there is an important difference in the processing of the coarsest scale. Since the coarsest scale of the NILC filter is not localized, a plain NILC map would be equivalent to a pixelbased ILC for all the multipoles of that scale. This procedure, however, is known to be quite susceptible to the “ILC bias”, due to chance correlations between the CMB and foregrounds. In order to mitigate this effect, the (single) covariance matrix which determines the ILC coefficients at the coarsest scale is not computed as a pixel average, but is rather estimated in the harmonic domain as an average over spherical harmonic coefficients using a spectral weight which equalizes the power of the CMB modes (based on a fiducial spectrum). This can be shown significantly to decrease the large scale errors.
In practice, our NILC processing depends on several implementation choices, as follows:

Input channels: in this work, the NILC algorithm isapplied to all Planck channelsfrom 44 to 857 GHz, omitting onlythe 30 GHz channel.

Preprocessing of point sources: identical to the SMICA preprocessing (see Appendix D).

Masking and inpainting: the NILC CMB map is actually produced in a threestep process. In a first step, the NILC weights are computed from covariance matrices evaluated using a Galactic mask removing about 2% of the sky (and is apodized at 1°). In a second step, those NILC weights are applied to needlet coefficients computed over the complete sky (except for point source masking/subtraction), yielding a NILC CMB estimate over the full sky (except for the point source mask). In short, the weights are computed over a masked sky but are applied to a full sky (up to point sources). In a final step, the pixels masked due to point source processing are replaced by the values of a constrained Gaussian realization (inpainting).

Spatial localization: the boundaries of the zones used for spatial localisation (shown at Fig. B.1) are obtained as isolevel curves of a low resolution map of Galactic emission.

Beam and transfer function: as in the SMICA processing, the resolution of the input maps is adjusted so that they have a beam of 5′. Therefore the resulting CMB map is automatically synthesized with an effective Gaussian beam of 5′, according to the unbiased nature of the ILC.

Using SMICA recalibration: in our current implementation, the NILC solution uses the values determined by SMICA for the CMB spectrum, given in Eq. (D.6).
Appendix C: Template fitting
The original SEVEM algorithm produced clean CMB maps at several frequencies through template fitting, followed by an estimation of the CMB power spectrum from these clean maps using a method based on the Expectation Maximization algorithm (MartınezGonzalez et al. 2003; Leach et al. 2008; FernándezCobos et al. 2012). From this power spectrum, a multifrequency Wienerfiltered CMB map was produced. For the present work, only the first step of the method, producing clean CMB maps at different frequencies, is considered. In addition, two of these clean maps are optimally combined to produce a final CMB map.
The templates used for cleaning are internal, i.e., they are constructed from Planck data, avoiding the need for external data sets, which usually complicate the analyses and may introduce inconsistencies. In the cleaning process, no assumptions about the foregrounds or noise levels are needed, rendering the technique very robust. The fitting can be done in real or wavelet space (using a fast wavelet adapted to the HEALPix pixelization; Casaponsa et al. 2011) to properly deal with incomplete sky coverage. For expediency, however, we fill in the small number of unobserved pixels at each channel with the mean value of their neighbouring pixels before applying SEVEM.
Linear coefficients, α_{j}, and templates used to clean individual frequency maps with SEVEM.
We construct our templates by subtracting two close Planck frequency channel maps, after first smoothing them to a common resolution, and converting to CMB temperature units if necessary, to ensure that the CMB signal is properly removed. A linear combination of the templates, t_{j}, is then subtracted from (hitherto unused) map, d, to produce a clean CMB map at that frequency. This is done either in real or wavelet space (i.e., scale by scale) at each position on the sky, (C.1)where n_{t} is the number of templates. If the cleaning is performed in real space, the α_{j} coefficients are obtained by minimizing the variance of the clean map, T_{c}, outside a given mask. When working in wavelet space, the cleaning is performed in the same way at each wavelet scale independently (i.e., the linear coefficients depend on the scale). Although we exclude very contaminated regions during the minimization, the subtraction is performed for all pixels and, therefore, the cleaned maps cover the fullsky (although foreground residuals are expected to be present in the excluded areas).
An additional level of flexibility may also be considered: the linear coefficients can be fixed over the full sky, or in several regions. The regions are then combined in a smooth way, by weighting the pixels at the boundaries to reduce possible discontinuities in the clean maps.
Since the method is linear, we may easily propagate the noise properties to the final CMB map. Moreover, it is very fast and permits the generation of thousands of simulations to characterize the statistical properties of the outputs, a critical need for many cosmological applications. The final CMB map retains the angular resolution of the original frequency map.
There are several configurations possible for SEVEM, depending the number of frequency maps to be cleaned or the number of templates used in the fitting. Note that the production of clean maps at different frequencies is of great interest in order to test the robustness of the results. Therefore, to define the best strategy, one needs to find a compromise between the number of maps that can be cleaned independently and the number of templates that can be constructed.
In particular, we have cleaned the 143 GHz and 217 GHz maps using four templates constructed as the difference of the following Planck channels (smoothed to a common resolution): (30−44), (44−70), (545−353), and (857−545). For simplicity, the two maps have been cleaned in real space, since there was no significant improvement found when using wavelets, especially at high latitude. In order to take into account the different spectral behaviour of the foregrounds at low and high Galactic latitudes, we considered two independent sky regions, using different sets of coefficients (see Table C.1 for the values of the linear coefficients for the main considered region). The first region corresponds to the brightest 3% of Galactic emission, while the second region is defined by the remaining 97% of the sky (see Sect. 4 for a detailed description of our Galactic mask construction). For the first region, the coefficients are actually estimated over the complete sky (we find that this is better than performing the minimization on only the brightest 3% of the sky, where the CMB is very subdominant), while for the second region, we exclude the bright 3% sky fraction, point sources detected at any frequency and those pixels which have not been observed in all channels.
Note that, for consistency, we have used the same configuration (four templates, cleaning in real space, two regions) for the analysis of the FFP6 simulations. However, we find that this simple configuration produces more contaminated CMB maps than for the data (although the region outside the confidence mask still has low contamination), indicating some differences between the foreground level in the data and in simulations. Therefore, conclusions derived from the FFP6 results for SEVEM should be taken with caution, since they are expected to provide overestimated residuals.
Our final CMB map was constructed by combining the 143 and 217 GHz maps by weighting the maps in harmonic space, taking into account the noise level, the resolution, and a rough estimation of the foreground residuals of each map (obtained from realistic simulations). This final map has a resolution corresponding to a Gaussian beam of 5′ FWHM.
Moreover, additional clean CMB maps (at frequencies 44, 70, 100, and 353 GHz) were also produced using different combinations of templates. In particular, to clean the 100 GHz map, we used the same templates and regions as for 143 and 217 GHz. This allowed us to produce three (almost) independent clean maps, in the sense that none of the three maps to be cleaned was used to construct the templates. For 44, 70, and 353 GHz, different combinations of templates were used, and the linear coefficients were chosen to be the same over the full sky. They were obtained by minimizing the variance of the map outside the same mask as that used to clean the central frequency maps on the largest region. The templates and the corresponding linear coefficients used for each of the considered frequencies are given in Table C.1.
The SEVEM clean frequency maps have been used in analyses of the isotropy and statistics of the CMB (Planck Collaboration XXIII 2014) and to obtain cosmological constraints from the ISW effect (Planck Collaboration XIX 2014). In particular, clean maps from 44 to 353 GHz were used for the stacking analysis presented in Planck Collaboration XIX (2014), while frequencies from 70 to 217 GHz were used for consistency tests in Planck Collaboration XXIII (2014).
The confidence mask provided for SEVEM is constructed by combining four types of selected regions. In particular, it excludes zones with high residuals identified through the subtraction of SEVEM clean maps at different frequencies, as well as the sources detected at all Planck channels using the Mexican hat wavelet algorithm (Planck Collaboration XXVIII 2014). These point sources are masked with holes of varying radius, according to the flux of each source. The pixels that are not observed by all channels are also masked. Finally, to ensure that the area left outside the mask is statistically robust, we also exclude from the analysis the brightest 20% of Galactic emission, leaving a useful area of around 76%. This provides a conservative mask for CMB analysis; however, we point out that smaller masks could also be used in specific applications, such as the lensing potential reconstruction described in Sect. 7).
Appendix D: Spectral matching
SMICA (Spectral Matching Independent Component Analysis) reconstructs a CMB map as a linear combination in harmonic space of N_{chan} input frequency maps with weights that depend on multipole ℓ. Given the N_{chan} × 1 vector x_{ℓm} of spherical harmonic coefficients for the input maps, it computes coefficients ŝ_{ℓm} for the CMB map as (D.1)where the N_{chan} × 1 vector w_{ℓ} containing the multipoledependent weights is chosen to give unit gain to the CMB with minimum variance. This is achieved with (D.2)where vector a is the spectrum of the CMB evaluated at each channel (allowing for possible interchannel recalibration factors) and C_{ℓ} is the N_{chan} × N_{chan} spectral covariance matrix of x_{ℓm}. Taking C_{ℓ} in Eq. (D.2) to be the sample spectral covariance matrix of the observations, (D.3)would implement a simple harmonicdomain ILC, similar to Tegmark et al. (2003). At the largest scales, we instead use a model, C_{ℓ}(θ), and determine the covariance matrix to be used in Eq. (D.2) by fitting C_{ℓ}(θ) to . This is done in the maximum likelihood sense for stationary Gaussian fields, that is, the best fit matrices, , are obtained for (D.4)Equations (D.1)–(D.4) summarize the basic principles of SMICA; its actual operation depends on a choice for the spectral model C_{ℓ}(θ), and on several implementationspecific details, which we briefly describe below.
We model the data as a superposition of CMB, noise and foregrounds. The latter are not parametrically modelled; instead, we represent the total foreground emission by d templates with arbitrary frequency spectra, angular spectra and correlations. In the spectral domain, this is equivalent to modelling the covariance matrices as (D.5)where C_{ℓ} is the angular power spectrum of the CMB, A is a N_{chan} × d matrix, P_{ℓ} is a positive definite d × d matrix, and N_{ℓ} is a diagonal matrix representing the noise power spectra of the data. The parameter vector θ contains all or part of the quantities in (D.5).
The decomposition D.5 reflects the fact that CMB, foregrounds and noise are independent components of the signal. Thus, SMICA is an ICA (independent component analysis) method. It operates by matching the observations to the spectral model (D.5) using the criterion (D.4).
The maximal flexibility in a SMICA fit of model (D.5) is obtained with all the parameters free, that is without any constraint on the spectrum C_{ℓ}, on the diagonal entries of N_{ℓ}, on a, or on A and P_{ℓ}. One would ideally fit all those parameters (except for obvious degeneracies, like that between a scale factor in a and the overall normalization of the CMB spectrum C_{ℓ}) over the whole multipole range. In practice, this turns out to be too difficult, given the large dynamic range both over the sky and over multipoles. We resort to a pragmatic threestep approach in which the criterion (D.4) is minimized by first fitting a, then A, and finally the linear parameters C_{ℓ} and N_{ℓ}. Each fit is conducted over the multipole ranges and the sky fraction most appropriate for the parameter of interest, as follows.
We first estimate the CMB spectral law a by fitting all model parameters (that is, without constraint) over a clean fraction of sky (f_{sky} = 40 %) in the range 100 ≤ ℓ ≤ 680 where the signal is CMBdominated in most of the channels and the beam window functions are accurately known. In this fit, which is done over a clean part of the sky, we use a foreground emission matrix, A, with only four columns. From this step, we only retain the best fit value for vector a. In the second step, we estimate the foreground emissivity by fixing a to its value from the previous step and fitting all the other parameters over a large fraction of sky (f_{sky} = 97%) in the range 4 ≤ ℓ ≤ 150 where the signal is dominated by the Galactic emission in all channels. From this second step, we retain the best fit value for the matrix A which, again, is adjusted without constraint other than having d = 6 columns. In the last step, we fit all power spectrum parameters: we fix a and A to their previously found values and fit for C_{ℓ} and P_{ℓ} at each ℓ.
Note that the first step (fitting a) amounts to recalibrating the input maps on the basis of CMB anisotropies. For the maps in thermodynamics units, we find (D.6)The value at 857 GHz is not accurately recovered by SMICA, so we have set a_{857} = 1. Since the norm of a is degenerate with a global scale factor for the CMB angular spectrum, it can only be recovered by SMICA up to a scale factor. This degeneracy is fixed here by taking a_{143} = 1. The recalibration step could have been omitted, since is very close the unit vector. However, we found that using improved the behavior of SMICA over using a = [ 1,...,1 ].
Before describing implementation details, we explain how SMICA deals with the varying resolution of the input channels, since the discussion thus far assumed that all input maps had the same resolution. Since SMICA works in the harmonic domain, it is a simple matter to account for the beam transfer function, b_{i}(ℓ), of the ith input map. The CMB sky multipoles s_{ℓm} contribute s_{ℓm}a_{i}b_{i}(ℓ)p_{i}(ℓ) to the harmonic coefficient of the ith map (where p_{i}(ℓ) is the pixel window function for the HEALPix map at ). Therefore, in order to produce a final CMB map at 5′ resolution, close to the highest resolution of Planck, we only need to work with input spherical harmonics rebeamed to 5′; that is, to apply SMICA on vectors with entries , where b_{5}(ℓ) is a 5′ Gaussian beam function. By construction, SMICA then produces an CMB map with an effective Gaussian beam of 5′ (without the pixel window function).
We now give further details on the actual implementation of SMICA:

Inputs: SMICA uses all nine Planck frequency channels from 30 to 857 GHz, harmonically transformed up to ℓ = 4000.

Preprocessing of point sources: SMICA is applied on input maps in which point sources are subtracted or masked. We start by fitting the PCCS point sources with S/N> 5 to a Gaussian shape where the source amplitude is estimated, together with its position and a constant factor representing the background variance. If the fit is successful (χ^{2} ≤ 2), the fitted point source is removed from the map; otherwise it is masked in all channels and the hole is inpainted by a simple diffusive filling process. This is done at all frequencies except 545 and 857 GHz, where all point sources with S/N> 7.5 are masked and inpainted.

Beams: when the harmonic coefficients of the input maps are rebeamed at 5′, we do not apply exactly the expression mentioned above, because the factor 1 /b_{i}(ℓ) would diverge at high ℓ for the lowest resolution input channels. That may not be a problem in infinite precision arithmetic, but would lead to matrices with extremely large condition numbers. Instead, we rebeam with the factor 1 /b_{i}(ℓ) replaced by min(1 /b_{i}(ℓ),1000). The rebeaming of the CMB modes then is no longer perfect, but this is of course irrelevant because the thresholding occurs in a regime where the signal is completely dominated by the noise, so that the contribution of the corresponding channel is already highly attenuated by the SMICA weights (as shown in Fig. D.1).
Fig. D.1 Weights, w_{ℓ}, given by SMICA to the input maps, after they are rebeamed to 5′and expressed in K_{RJ}, as a function of multipole. Top panel: linear scale; bottom panel: the absolute value of the weights on a logarithmic scale.
Fig. D.2 Contribution of each input channel to the noise in the SMICA map.

Masking: in practice, SMICA operates on a masked sky, the mask being applied after the point source processing. The mask is obtained by thresholding a heavily smoothed version of the point source mask. The threshold is chosen to leave about 97% of the sky. Because of the heavy smoothing, the mask has smooth contours and is only sensitive to large aggregates of point sources: the masked areas mostly lie in the Galactic plane, but include also a few bright regions like the Large Magellanic Cloud.

Inpainting: the SMICA map used in this paper has no real power in the masked region described above. However, for convenience, an inpainted SMICA map has also been produced by replacing the masked pixels with a constrained Gaussian realization obtained by the method of BenoitLévy et al. (2013). That map appears in Planck Collaboration I (2014).

Binning: in our implementation, we use binned spectra.

Processing at fine scales: since there is little point trying to model the spectral covariance at high multipoles, because the sample estimate is sufficient, SMICA implements a simple harmonic ILC at ℓ > 1500; that is, it applies the filter (D.2) with .

Confidence mask: a confidence mask (Fig. 3) is provided with SMICA, constructed in the following way. The SMICA CMB map is bandpass filtered through a spectral window v(ℓ) = exp [−((ℓ − 1700)/200)^{2}/2 ]. The result is squared and smoothed at 2° resolution, yielding a map of the (bandpassed) variance of the CMB map. That variance is corrected for the noise contribution by subtracting the variance map for the noise obtained by the same procedure applied to the SMICA HRHD map. If the SMICA map contained only CMB and noise, the variance map would have a uniform value ∑ _{ℓ}v(ℓ)^{2}b_{5}(ℓ)^{2}C(ℓ)(2ℓ + 1)/4π = 31.3 μK^{2} over the sky. The confidence map is obtained by thresholding the noisecorrected variance map at 70 μK^{2}.
Viewed as a filter, SMICA can be summarized by the weights w_{ℓ} applied to each input map as a function of multipole. In this sense, SMICA is strictly equivalent to coadding the input maps after convolution by specific axisymmetric kernels directly related to the corresponding entry of w_{ℓ}.
The SMICA weights used here are shown in Fig. D.1 (for input maps in units of K_{RJ}). We see, in particular, the (expected) progressive attenuation of the lowest resolution channels with increasing multipole. Figure D.2 shows the contribution of each input channel to the noise in the SMICA CMB map as a function of multipole. The spectral noise contribution from channel i is simply obtained as w(ℓ)^{2}N_{i}(ℓ), where w_{i}(ℓ) is the ith entry of the weight vector w(ℓ) and N_{i}(ℓ) is the angular spectrum of the ith noise map.
More details about SMICA are given in Cardoso et al. (2008), as well as in applications to the analysis of WMAP (Patanchon et al. 2005) and Archeops data (Tristram et al. 2005). An application to the measurement of the tensortoscalar ratio using CMB Bmodes is discussed in Betoule et al. (2009). Within the Planck collaboration, SMICA is used to define the Plik highℓ likelihood (Planck Collaboration XV 2014), but physical models of foreground emission are used there instead of the nonparametric foreground model used here. SMICA is also used to crosscheck the HFI calibration (see Planck Collaboration VIII 2014).
Appendix E: FFP6 simulation results
Fig. E.1
Standard deviation of the four foregroundcleaned CMB maps from the FFP6 simulation. All maps have been downgraded to a HEALPix resolution of N_{side} = 128. 
Fig. E.2
Pairwise differences between foregroundcleaned CMB maps from the FFP6 simulation. All maps have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale differences. 
Fig. E.3
Angular power spectra of residual foreground emission in the CMB maps from the FFP6 simulation. The components shown are: thermal dust; cosmic infrared background fluctuations (“firb”); point sources; CMB anisotropies; noise; and the sum of all others. From top to bottom, the panels show the results for CommanderRuler, NILC, SEVEM, and SMICA. 
Here we study the performance of the component separation algorithms on the FFP6 simulation described in Sect. 4, providing additional information beyond that in the body of the paper. Much of the analysis presented here mirrors that shown for the data in Sects. 5 and 6, allowing a direct comparison between the two analyses.
Appendix E.1: CMB maps
First, we show in Fig. E.1 the standard deviation between the four foregroundcleaned FFP6 maps, similar to that shown in Fig. 5 for the data. Figure E.2 shows all pairwise differences between the same maps, mirroring Fig. 6 for the data. These two plots highlight an important point concerning the FFP6 analysis already mentioned in Sect. 5, namely that in nearGalactic regions, where the foregrounds are important, the internal differences between the four algorithms are larger in the FFP6 simulation than in the real data. This is due to the fact that each component separation algorithm has been optimized using the real data in terms of model definition, localization, etc. Then, the same models have been used for the FFP6 simulation without change. Only the parameters within those models are refitted to the new data set. This implies, in fact, that we expect each method to perform better on the data than the simulations in terms of absolute residuals, to the extent that the simulation matches the real sky. In other words, the FFP6 simulation provides a conservative estimate of the residual errors in the real data.
Appendix E.2: Power spectrum residuals from individual components
In Fig. E.3 we show the residual effect of some of the individual components on the foregroundcleaned CMB map. The thermal dust emission, CIB fluctuations, point sources, and noise have been processed individually with each algorithm. All other components (freefree, synchrotron, spinning dust, CO, thermal SZ, and kinetic SZ) are shown as a single, composite residual component.
Appendix E.3: CMB power spectra and cosmological parameters
Fig. E.4
Estimates of the CMB power spectra from the foregroundcleaned FFP6 maps, computed by XFaster. The solid lines show the spectra after subtracting the bestfit model of residual foregrounds. The vertical dotted line shows the maximum multipole (ℓ = 2000) used in the likelihood for fitting the foreground model and cosmological parameters (see Sect. E.3 for further details). The dashed lines show the spectra before residual foreground subtraction. 
We assess the performance of our component separation techniques by evaluating cosmological constraints from the foregroundcleaned CMB maps derived from the FFP6 simulation.
Figure E.4 shows the estimates of the angular power spectra of the CMB maps. Figure E.5 compares the cosmological parameters derived from the four foregroundcleaned CMB maps, together with CamSpec^{7} and Plik, to the input (theoretical) parameters for different ℓranges. The parameter space is defined by the same model applied to the real data in Sect. 6, including six ΛCDM and two foreground parameters. All deviations from input parameters are small and within 1σ up to ℓ = 2000, verifying that all methods work well in this multipole range. However, for ℓ_{max} = 2500 we start to see significant shifts, e.g., for Ω_{b}h^{2} and n_{s}. Further, the point source foreground parameter, A_{PS}, reaches large values, implying that assumptions concerning the highℓ foreground model become important. For these reasons, we consider ℓ_{max} = 2000 as the maximum recommended ℓ range for these maps in the current data release.
Still, the overall agreement is excellent between all codes and all ℓ ranges. In particular, we see that differences in the bandpower spectra at high ℓ between the different codes are mostly absorbed by the twoparameter foreground model. For instance, the CommanderRuler band power spectrum has more power at high ℓ due to noise or residual point sources, but this excess is well fitted by the twoparameter foreground model, and mostly interpreted in terms of a residual point source component; this is expected, given the lower angular resolution of this map.
As mentioned above, n_{s} and A_{s} are to some extent sensitive to ℓ_{max}. These parameters are degenerate with the foreground parameters. This may suggest that our C_{ℓ} foreground templates deviate more from the shape of the Poissonian and clustered component in the CMB map. This is a limitation of the simple foreground templates used here. To properly describe the foreground residuals in the reconstructed maps, we should use a foreground power spectrum template tailored to each method. For instance, such templates may be constructed by processing simulated foreground maps though each of the four pipelines. The templates are then given by the pseudoC_{ℓ} of each of the processed foreground map. However, our analysis shows that the current simple model provides accurate results when restricting the analysis to ℓ_{max} = 2000.
Figure E.6 shows the bestfit power spectrum residuals for the CMB map, CamSpec and Plik relative to the input CMB ΛCDM model estimated up to ℓ = 2000. These plots show that the residuals of the CMB mapbased bestfit models are comparable to the CamSpec and Plik residuals, and smaller than 40 μK^{2} for most of the ℓ range, with larger deviations observed for CamSpec at ℓ ~ 200. At higher ℓs the residuals are smaller than 10 μK^{2} for both approaches, all showing similar trends. Thus, both the map and spectrumbased likelihoods recover input parameters reasonably well, with the latter yielding slightly larger deviations from the bestfit model of the input CMB realization.
Fig. E.5
Comparison of cosmological parameters derived from the FFP6 simulation using different methods. The parameters shown as blue, red and green points indicate results obtained with ℓ_{max} = 1500, 2000 and 2500, respectively, and the yellow points show the results derived by CamSpec and Plik using crossspectra. The black horizontal lines mark the input parameter values. The values of the foreground parameters are not shown for CamSpec or Plik since they use a different model. The matter power spectrum pivot scale was k = 0.002 for all likelihoods, except CamSpec for which k = 0.05 was used. 
Fig. E.6
Residuals of mapbased and spectrumbased bestfit models relative to the FFP6 simulation input ΛCDM spectrum, shown for each algorithm up to ℓ_{max} = 2000. Cosmic variance is shown as the black dashed line. 
All Tables
Estimated monopoles and dipoles in Galactic coordinates, all measured in thermodynamic μK.
Linear coefficients, α_{j}, and templates used to clean individual frequency maps with SEVEM.
All Figures
Fig. 1
Foregroundcleaned CMB maps derived by CommanderRuler, NILC, SEVEM and SMICA. Note that the SMICA map has been filled in smoothly inside a 3% Galactic mask. 

In the text 
Fig. 2
Combined Galactic (CG) emission masks for the Planck data, corresponding to sky fractions of 20, 40, 60, 70, 75, 80, 90, 97, and 99%. The masks are named CG20, etc. 

In the text 
Fig. 3
Summary of component separation (CS) confidence masks. Each pixel is encoded in terms of a sum in which CommanderRuler equals 1 (light blue), NILC equals 2 (dark red), SEVEM equals 4 (yellow), and SMICA equals 8 (light red). The masks are named CSCR75, CSNILC93, CSSEVEM76, and CSSMICA89, respectively, reflecting their accepted sky fraction. The union mask (U73), used for evaluation purposes in this paper, removes all coloured pixels. 

In the text 
Fig. 4
Beam transfer functions of the four foregroundcleaned CMB maps. 

In the text 
Fig. 5
Standard deviation between the four foregroundcleaned CMB maps. All maps have been downgraded to a HEALPix resolution of N_{side} = 128. The differences are typically less than 5 μK at high Galactic latitudes, demonstrating that the maps are consistent over a large part of the sky. 

In the text 
Fig. 6
Pairwise differences between foregroundcleaned CMB maps. All maps have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale differences. The linelike discontinuities in the differences involving SEVEM is due to the two different regions used in this algorithm to clean the sky (see Appendix C for details). 

In the text 
Fig. 7
CMB residual maps from the FFP6 simulation. A monopole determined at high Galactic latitude has been subtracted from the maps, and they have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale features. The residuals presented here provide a conservative estimate of those expected in the data (see text for details). 

In the text 
Fig. 8
Angular power spectra of the foregroundcleaned CMB maps and halfring halfdifference (HRHD) maps. The spectra have been evaluated using the U73 mask apodized with a 30′ cosine function. 

In the text 
Fig. 9
Angular power spectra of FFP6 simulated components evaluated over the common mask (top) and the common point source mask (bottom), both apodized with a 30′ cosine function. Three components are shown: the CMB (dashed line); noise (dotdashed line); and the sum of all foregrounds (solid line). A nonlinear scale is used on the horizontal axis to show all the features of the spectra. 

In the text 
Fig. 10
Estimates of the CMB power spectra from the foregroundcleaned maps, computed by XFaster. The solid lines show the spectra after subtracting the bestfit model of residual foregrounds. The vertical dotted line shows the maximum multipole (ℓ = 2000) used in the likelihood for fitting the foreground model and cosmological parameters (see Sect. 6.2.2 for further details). The dashed lines show the spectra before residual foreground subtraction. 

In the text 
Fig. 11
Comparison of cosmological and foreground parameter values estimated from the foregroundcleaned CMB maps for ℓ_{max} = 2000 (in red) and those obtained with CamSpec and Plik likelihoods (in blue). The values of the foreground parameters are not shown for CamSpec and Plik, since they use a different foreground model. 

In the text 
Fig. 12
Residuals of all mapbased bestfit models relative to CamSpec bestfit model (assuming a prior on τ) for ℓ_{max} = 2000. 

In the text 
Fig. 13
Lowℓ power spectrum amplitude and tilt constraints measured relative to the bestfit PlanckΛCDM model derived from foregroundcleaned CMB maps smoothed to 6° FWHM. The cross shows the bestfit model (q,n) = (1,0). 

In the text 
Fig. 14
Lensing power spectrum estimates from FFP6 simulations using an apodized mask covering f_{sky,2} ≃ 0.70 of the sky. 

In the text 
Fig. 15
Posterior mean foreground amplitude maps derived from the lowresolution analysis. From top to bottom are shown the lowfrequency, CO and thermal dust emission maps. 

In the text 
Fig. 16
Posterior mean spectral parameter maps derived from the lowresolution analysis. The top panel shows the power law index of the lowfrequency component, and the bottom panel shows the emissivity index of the onecomponent thermal dust model. Note that the systematic error due to monopole and dipole uncertainties is significant for the dust emissivity in regions with a low thermal dust amplitude. 

In the text 
Fig. 17
χ^{2} per pixel for the joint CMB and foreground analysis. The expected value for an acceptable fit is 7, corresponding to the number of frequency bands used in this analysis. The pixels with high values can be classified into two types, due to either modelling errors (i.e., high residuals in the Galactic plane) or to unmodelled correlated noise (i.e., stripes crossing through low dust emission regions). 

In the text 
Fig. 18
Comparison of the highresolution Ruler (top) and lowresolution Commander (bottom) amplitude maps for a particularly strong CO complex near the Fan region; the maps are centred on Galactic coordinates (l,b) = (110°,15°), and the grid spacing is 5°. Columns show, left to right: the CMB amplitude; the CO amplitude at 100 GHzand the thermal dust amplitude at 353 GHz. 

In the text 
Fig. 19
Amplitude residual maps, A_{out} − A_{in}, computed blindly from the FFP6 simulation. The panels show (from top to bottom) the lowfrequency residual at 30 GHz, the CO residual at 100 GHz and the thermal dust residual at 353 GHz. All units are thermodynamic μK. The white lines indicate the boundary of the Commander likelihood analysis mask, removing 13% of the sky. 

In the text 
Fig. 20
Error validation for component amplitudes, evaluated from the FFP6 simulation. The upper panel shows histograms of the normalized errors δ = (A_{out} − A_{in}) /σ_{out} for the three foreground components and including all pixels outside the Commander likelihood analysis mask. The lower panel shows histograms of the fractional error f ≡ (A_{out} − A_{in}) /A_{in} for pixels with a foreground detection level above 5σ. No evidence of significant bias is observed for any component, and the uncertainty estimates for the lowfrequency and CO components are accurate to about 12%; the thermal dust uncertainty is underestimated by a factor of 2 due to the presence of unmodelled fluctuations. 

In the text 
Fig. 21
Validation of spectral parameters for lowfrequency foregrounds, thermal dust, and CO emission, evaluated from the FFP6 simulation. Each histogram shows the error distribution at the two leading subdominant frequencies in the form of the normalized errors δ = (A_{out}(ν) − A_{in}(ν)) /σ_{out}(ν) for all pixels outside the Commander likelihood analysis mask, where A_{out}(ν) is the predicted foreground amplitude at frequency ν given the estimated amplitude and spectral parameters, and σ_{out}(ν) is the corresponding standard deviation computed over the sample set. 

In the text 
Fig. 22
Top: difference map between the thermal dust amplitude at 353 GHz presented in this paper using the Planck 30 to 353 GHz frequencies, and that derived by Planck Collaboration XI (2014) using only the Planck 353, 545 and 857 GHz channels and the 100 μm IRIS map. Both maps are smoothed to a common resolution of 40′. Middle: the same difference map, but accounting for relative monopole, dipole and zodiacal emission treatment. Bottom: T–T plot between the same two maps after applying relative corrections. Any pixels with a CO amplitude at 100 GHz larger than 1 K km s^{1} are removed from this plot. 

In the text 
Fig. 23
Top: dominant foreground component per pixel at 30 GHz in the FFP6 simulation. Dark blue indicates that synchrotron emission is the strongest component at 30 GHz, light blue indicates that freefree dominates, and orange indicates that spinning dust (AME) is the strongest component. Bottom: the recovered lowfrequency powerlaw index derived from the same simulation. 

In the text 
Fig. B.1
Spectral localization for NILC using nine spectral window functions defining nine needlet scales (top panel). The scaledependent spatial localization partitions the sky into one zone (for scale 1), two zones (for scale 2), four zones (for scale 3), or twenty zones (for scales 5, 6, 7, 8 and 9). The twozone, 4zone and 20zone partitions are depicted in the lower panels. 

In the text 
Fig. D.1
Weights, w_{ℓ}, given by SMICA to the input maps, after they are rebeamed to 5′and expressed in K_{RJ}, as a function of multipole. Top panel: linear scale; bottom panel: the absolute value of the weights on a logarithmic scale. 

In the text 
Fig. D.2
Contribution of each input channel to the noise in the SMICA map. 

In the text 
Fig. E.1
Standard deviation of the four foregroundcleaned CMB maps from the FFP6 simulation. All maps have been downgraded to a HEALPix resolution of N_{side} = 128. 

In the text 
Fig. E.2
Pairwise differences between foregroundcleaned CMB maps from the FFP6 simulation. All maps have been downgraded to a HEALPix resolution of N_{side} = 128 to show the largescale differences. 

In the text 
Fig. E.3
Angular power spectra of residual foreground emission in the CMB maps from the FFP6 simulation. The components shown are: thermal dust; cosmic infrared background fluctuations (“firb”); point sources; CMB anisotropies; noise; and the sum of all others. From top to bottom, the panels show the results for CommanderRuler, NILC, SEVEM, and SMICA. 

In the text 
Fig. E.4
Estimates of the CMB power spectra from the foregroundcleaned FFP6 maps, computed by XFaster. The solid lines show the spectra after subtracting the bestfit model of residual foregrounds. The vertical dotted line shows the maximum multipole (ℓ = 2000) used in the likelihood for fitting the foreground model and cosmological parameters (see Sect. E.3 for further details). The dashed lines show the spectra before residual foreground subtraction. 

In the text 
Fig. E.5
Comparison of cosmological parameters derived from the FFP6 simulation using different methods. The parameters shown as blue, red and green points indicate results obtained with ℓ_{max} = 1500, 2000 and 2500, respectively, and the yellow points show the results derived by CamSpec and Plik using crossspectra. The black horizontal lines mark the input parameter values. The values of the foreground parameters are not shown for CamSpec or Plik since they use a different model. The matter power spectrum pivot scale was k = 0.002 for all likelihoods, except CamSpec for which k = 0.05 was used. 

In the text 
Fig. E.6
Residuals of mapbased and spectrumbased bestfit models relative to the FFP6 simulation input ΛCDM spectrum, shown for each algorithm up to ℓ_{max} = 2000. Cosmic variance is shown as the black dashed line. 

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.