Issue 
A&A
Volume 596, December 2016



Article Number  A109  
Number of page(s)  26  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201629022  
Published online  12 December 2016 
Planck intermediate results
XLVIII. Disentangling Galactic dust emission and cosmic infrared background anisotropies
^{1} APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de ParisSorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet 75205 Paris Cedex 13 France
^{2} African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg 7945, Cape Town, South Africa
^{3} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{4} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{5} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{6} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{7} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{8} California Institute of Technology, Pasadena, California, CA 91125, USA
^{9} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, CA 94720, USA
^{10} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{11} Département de Physique Théorique, Université de Genève, 24 Quai E. Ansermet, 1211 Genève 4, Switzerland
^{12} Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{13} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33007 Oviedo, Spain
^{14} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{15} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{16} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles, CA 90089, USA
^{17} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{18} Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{19} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00014 Helsinki, Finland
^{20} Department of Physics, Princeton University, Princeton, New Jersey, NJ 08544, USA
^{21} Department of Physics, University of California, Santa Barbara, California, CA 93106, USA
^{22} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{23} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{24} Dipartimento di Fisica e Astronomia, Alma Mater Studiorum, Università degli Studi di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{25} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{26} Dipartimento di Fisica, Università La Sapienza, P. le A. Moro 2, 00185 Roma, Italy
^{27} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy
^{28} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, 34127 Trieste, Italy
^{29} Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, Roma, Italy
^{30} Dipartimento di Matematica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, 00133 Roma, Italy
^{31} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain
^{32} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{33} Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’ Aquila, Italy
^{34} HGSFP and University of Heidelberg, Theoretical Physics Department, Philosophenweg 16, 69120 Heidelberg, Germany
^{35} Helsinki Instituteof Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00014 Helsinki, Finland
^{36} INAF–Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{37} INAF–Osservatorio Astronomico di Roma, via di Frascati 33, 00040 Monte Porzio Catone, Italy
^{38} INAF–Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 40127 Trieste, Italy
^{39} INAF/IASF Bologna, via Gobetti 101, 40129 Bologna, Italy
^{40} INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy
^{41} INFN – CNAF, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{42} INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{43} INFN, Sezione di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{44} INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{45} INFN, Sezione di Roma 2, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, 00185 Roma, Italy
^{46} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{47} Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{48} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98 bis Boulevard Arago, 75014 Paris, France
^{49} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{50} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0371 Oslo, Norway
^{51} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, 38205 Tenerife, Spain
^{52} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain
^{53} Istituto Nazionale di Fisica Nucleare, Sezione di Padova, via Marzolo 8, 35131 Padova, Italy
^{54} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, CA 31109, USA
^{55} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{56} Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA
^{57} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{58} LAL, Université ParisSud, CNRS/IN2P3, 91898 Orsay, France
^{59} LERMA, CNRS, Observatoire de Paris, 61 Avenue de l’Observatoire, 75014 Paris, France
^{60} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{61} Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{62} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{63} Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{64} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{65} Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{66} Nicolaus Copernicus Astronomical Center, Bartycka 18, 00716 Warsaw, Poland
^{67} Nordita (Nordic Institute for Theoretical Physics), Roslagstullsbacken 23, 106 91 Stockholm, Sweden
^{68} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{69} School of Chemistry and Physics, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{70} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{71} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{72} Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby BC, Canada
^{73} Sorbonne UniversitéUPMC, UMR 7095, Institut d’Astrophysique de Paris, 98 bis Boulevard Arago, 75014 Paris, France
^{74} Space Sciences Laboratory, University of California, Berkeley, California, CA 94720, USA
^{75} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{76} The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{77} UPMC Univ Paris 06, UMR 7095, 98 bis Boulevard Arago, 75014 Paris, France
^{78} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{79} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{80} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
^{⋆}
Corresponding author: M. Remazeilles, email: mathieu.remazeilles@manchester.ac.uk
Received: 27 May 2016
Accepted: 9 August 2016
Using the Planck 2015 data release (PR2) temperature maps, we separate Galactic thermal dust emission from cosmic infrared background (CIB) anisotropies. For this purpose, we implement a specifically tailored componentseparation method, the socalled generalized needlet internal linear combination (GNILC) method, which uses spatial information (the angular powerspectra) to disentangle the Galactic dust emission and CIB anisotropies. We produce significantly improved allsky maps of Planck thermal dust emission, with reduced CIB contamination, at 353, 545, and 857 GHz. By reducing the CIB contamination of the thermal dust maps, we provide more accurate estimates of the local dust temperature and dust spectral index over the sky with reduced dispersion, especially at high Galactic latitudes above b = ±20°. We find that the dust temperature is T = (19.4 ± 1.3) K and the dust spectral index is β = 1.6 ± 0.1 averaged over the whole sky, while T = (19.4 ± 1.5) K and β = 1.6 ± 0.2 on 21% of the sky at high latitudes. Moreover, subtracting the new CIBremoved thermal dust maps from the CMBremoved Planck maps gives access to the CIB anisotropies over 60% of the sky at Galactic latitudes b > 20°. Because they are a significant improvement over previous Planck products, the GNILC maps are recommended for thermal dust science. The new CIB maps can be regarded as indirect tracers of the dark matter and they are recommended for exploring crosscorrelations with lensing and largescale structure optical surveys. The reconstructed GNILC thermal dust and CIB maps are delivered as Planck products.
Key words: cosmology: observations / methods: data analysis / ISM: general / dust, extinction / infrared: diffuse background / largescale structure of Universe
© ESO, 2016
1. Introduction
The various populations of dust grains in the Galaxy are heated by absorbing the ultraviolet emission from stars. By reemitting the light at infrared frequencies, the heated dust grains are responsible for the thermal dust radiation of the Galaxy. For this reason, the dust emission is a tracer of the gas and particle density in the interstellar medium (Planck Collaboration XIX 2011; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XXVIII 2015; Planck Collaboration Int. XXXI 2016) and of the star formation activity in the Galaxy (Draine & Li 2007). The Galactic thermal dust emission is also one of the major astrophysical foregrounds for observations of the cosmic microwave background (CMB) (Planck Collaboration Int. XXX 2016; BICEP2/Keck Array and Planck Collaborations 2015). Incorrect modelling of thermal dust might be responsible for a significant bias on the cosmological parameters (e.g., Remazeilles et al. 2016). The characterization of Galactic thermal dust emission over the whole sky is therefore essential for the accurate subtraction of this foreground from the CMB observations. Accurate characterization of the Galactic dust is also useful for the analysis of supernovae observations, where the Galactic dust causes extinction (Riess et al. 1996).
An unresolved background of dusty starforming early galaxies also generates diffuse emission, known as the cosmic infrared background radiation (Puget et al. 1996; Gispert et al. 2000; Lagache et al. 2005). The cosmic infrared background (CIB) anisotropies are a probe of star formation history in the Universe and also an indirect tracer of the dark matter (Planck Collaboration XVIII 2011; Planck Collaboration XXX 2014). At high frequencies (≳300 GHz), the Galactic thermal dust emission and the CIB radiation both scale approximately as modified blackbodies. This makes it challenging to separate the dust and CIB components solely on the basis of their spectral properties (Planck Collaboration XI 2014; Planck Collaboration X 2016)
The previously released Planck^{1} dust maps – the Planck 2013 (P13) dust model (Planck Collaboration XI 2014) and the Planck 2015 (P15) dust model (Planck Collaboration X 2016) – have been produced by fitting a modified blackbody (MBB) spectrum to the Planck data. For the P13 dust map, a standard χ^{2} fitting of the MBB spectrum was applied pixel by pixel to four maps, namely the CMBremoved Planck temperature maps at 353, 545, and 857 GHz from the Planck 2013 data release (hereafter PR1), and a 100 µ map obtained from a combination of the IRIS map from MivilleDeschênes & Lagache (2005) and the map from Schlegel et al. (1998). The CMB removal in the Planck frequency channels was performed by subtracting the PlanckSMICA CMB map (Planck Collaboration XII 2014) from the Planck frequency maps. For the P15 dust map, a Bayesian fitting of the MBB spectrum was implemented on the Planck 2015 data release (hereafter PR2) temperature maps by using the full set of Planck frequency channels.
However, the Planck dust models P13 and P15 still suffer from contamination by the CIB anisotropies. In particular, at high Galactic latitudes the contamination by CIB anisotropies adds significant uncertainty to the measured dust spectral index and dust temperature (Planck Collaboration XI 2014).
By definition, the spectral fitting employed in Planck Collaboration XI (2014) and Planck Collaboration X (2016) relied solely on the frequency information to reconstruct the Galactic thermal dust model from observations of the sky emission. Because the Galactic dust emission and the extragalactic CIB emission have such similar spectral indices in the Planck bands, the result of these frequencybased fits is that the CIB anisotropies inevitably leak into the Planck dust model maps. In order to disentangle the Galactic thermal dust emission and the extragalactic CIB emission, additional discriminating statistical information is required.
The CIB temperature fluctuations have been successfully measured by Planck in relatively small regions of the sky, where the Galactic dust contamination is low (Planck Collaboration XXX 2014), and the angular power spectra of the CIB anisotropies have been computed at frequencies from 143 GHz to 3000 GHz.
In this work we perform the separation of the Galactic thermal dust and CIB components over a large area of the sky by exploiting not only the frequency spectral information, but also the spatial information through the use of the Planck CIB bestfit angular power spectra computed in Planck Collaboration XXX (2014, hereafter CIB 2013). The CIB power spectrum scales approximately as ℓ^{1} (Planck Collaboration XVIII 2011), while the dust power spectrum scales approximately as ℓ^{2.7} (Planck Collaboration XXX 2014). This distinct spatial behaviour provides the necessary extra statistical information that enables robust separation of thermal dust emission and CIB radiation. Although the CIB 2013 angular power spectra have only been estimated in small areas of the sky, we assume that the statistics of the CIB anisotropies are the same in a larger area of the sky, because of the homogeneity and isotropy of the CIB emission.
The componentseparation method employed in this work is the generalized needlet internal linear combination (GNILC) method, first developed in Remazeilles et al. (2011b). It is worth noting that the GNILC method has also been applied in a different context to simulations of a radio intensity mapping experiment for separating the cosmological H i21cm temperature fluctuations and the Galactic synchrotron radiation in Olivari et al. (2016), where the componentseparation problem was similar.
This paper is organized as follows. In Sect. 2 we present the data used in the analysis. In Sect. 3 we give a summary of the componentseparation method that we implement on the data to disentangle the Galactic dust emission and the CIB anisotropies; the full description of the method and validation on simulations are presented in Appendix A. In Sect. 4 we discuss the results for the Galactic thermal dust emission and the estimated spectral parameters. In Sect. 5 we discuss the results for the CIB emission. In Sect. 6 we explore the correlations of the new dust and CIB maps with the H i gas tracer. We conclude in Sect. 7.
2. Data and preprocessing
2.1. Planck data
The data used in this paper are the temperature fullmission sky maps (Planck Collaboration VI 2016; Planck Collaboration VIII 2016) of the Planck 2015 data release (PR2) that have been made publicly available on the Planck Legacy Archive. We make use of the nine singlefrequency maps from 30 to 857 GHz from both LFI and HFI instruments. As discussed in Planck Collaboration XIV (2014), the zodiacal light emission is removed from the Planck HFI temperature maps (100 to 857 GHz) by fitting different Planck surveys of the sky with the COBE/DIRBE^{2} zodiacal model (Kelsall et al. 1998). Because different Planck surveys are taken at different times, the sky is observed through different column depths of interplanetary dust. Differencing two surveys removes all distant structure in the maps, such as Galactic and extragalactic emission, but leaves a detectable Zodiacal signal. This difference signal is fit to extend the COBE zodiacal model to Planck frequencies. The entire, undifferenced signal is then reconstructed from the model and removed from the data of each Planck HFI bolometer prior to mapmaking (Planck Collaboration VIII 2016).
We also make use of the Planck temperature halfmission sky maps (hereafter HM1 and HM2, as defined in Planck Collaboration VIII 2016) in the nine frequency channels, in order to estimate by their halfdifference (see Eq. (A.14)) the local rms of the instrumental noise in the Planck fullmission maps.
2.2. The IRAS 100 μm map
Following Planck Collaboration XI (2014), in addition to the Planck PR2 data we also use in this work the fullsky temperature map at 100 μm based on the combination of the IRIS map (MivilleDeschênes & Lagache 2005) and the map of Schlegel et al. (1998, hereafter SFD map) both projected on the HEALPix grid (Górski et al. 2005) at N_{side} = 2048. The combined 100 μm map is compatible with the SFD map at angular scales larger than 30′ and compatible with the IRIS map at smaller angular scales. The effective beam resolution of the combined 100 μm map is 4.3′ and the noise rms level is 0.06 MJy sr^{1}.
It should be noted that residual lowlevel zodiacal light emission is present in the combined 100 μm map because the zodiacal emission has not been corrected in the same way in the SFD map and in the IRIS map. We refer to the appendix of Planck Collaboration XI (2014) for further discussion of this point.
2.3. Preprocessing of the point sources
We make use of nine point source masks, one for each Planck frequency channel, in order to remove the point sources detected in each frequency at a signaltonoise ratio S/N> 5 in the second Planck Catalogue of Compact Sources, PCCS2 (Planck Collaboration XXVI 2016).
The masked pixels in each Planck frequency map are filled in through a minimum curvature spline surface inpainting technique, implemented in the Planck Sky Model (PSM) software package (Delabrouille et al. 2013) and described in Remazeilles et al. (2015) and Planck Collaboration XII (2016). For consistency, we also consider the sourcesubtracted version of the combined 100 μm map that is described in the appendix of Planck Collaboration XI (2014). The inpainted Planck 2015 maps and the inpainted combined 100 μm map are the inputs to the componentseparation algorithm described in the next section.
3. Summary of the componentseparation method
The componentseparation technique that we follow in this work is based on Remazeilles et al. (2011b) and called GNILC.
In order to simplify the reading of the paper, we give a brief summary of the method employed in this work to separate the thermal dust and CIB anisotropies. A complete description of the formalism and technical details of GNILC are presented in Appendix A.
Each frequency map is first decomposed on a needlet (spherical wavelet) frame (Narcowich et al. 2006; Guilloux et al. 2009). The localization properties of the needlets allow us to adapt the component separation to the local conditions of contamination in both harmonic space and real space (Delabrouille et al. 2009; Remazeilles et al. 2011b; Basak & Delabrouille 2012; Remazeilles et al. 2013; Basak & Delabrouille 2013). We define ten needlet windows, { h^{(j)}(ℓ) } _{1 ≤ j ≤ 10}, having a Gaussian shape and acting as bandpass filters in harmonic space, each of them selecting a specific subrange of angular scales (see Fig. A.2). The spherical harmonic transform, a_{ℓm}, of each frequency map is bandpass filtered in harmonic space by the ten needlet windows. The inverse transform of the bandpassfiltered coefficient, h^{(j)}(ℓ)a_{ℓm}, provides a needlet map at scale j, conserving only statistical information from the range of ℓ considered. Therefore, we have 100 input maps (10 frequencies times 10 needlet scales). The component separation is performed on each needlet scale independently. The main steps of the GNILC algorithm are the following. For each needlet scale, j, considered we perform seven steps:

1.
We compute the frequencyfrequency data covariance matrix,at pixel p, and scale j, (1)where is a domain of pixels centred at pixel p and and are the needlet maps at scale j of the observations for the pair of frequencies a,b. In practice, the domain of pixels, , is defined by the convolution in real space of the product of the needlet maps with a Gaussian kernel. The width of the Gaussian kernel is a function of the needlet scale considered.

2.
Similarly, we compute the frequencyfrequency covariance matrix of the instrumental noise, at pixel p, and scale j, (2)where the instrumental noise maps, n(p), are estimated from the halfdifference of the halfmission HM1 and HM2 Planck maps.

3.
Similarly, we compute the frequencyfrequency covariance matrix of the CMB, , and the frequencyfrequency covariance matrix of the CIB, , but this time using CMB maps and CIB maps that are simulated from the Planck CMB bestfit C_{ℓ} (Planck Collaboration XV 2014) and the Planck CIB bestfit (Planck Collaboration XXX 2014) respectively. The simulated maps were analyzed with the same needlet decomposition as was applied to the real data before computing the CMB and CIB covariance matrices.

4.
We compute the “nuisance” covariance matrix, , by coadding the noise covariance matrix, the CMB covariance matrix, and the CIB covariance matrix: (3)

5.
We diagonalize the transformed data covariance matrix (4)where N_{ch} is the number of frequency channels. In this representation, the eigenvalues that are close to unity correspond to the nuisance power (CIB plus CMB plus noise), while the m eigenvalues larger than unity that are collected in the diagonal matrix correspond to the power of the Galactic signal. The matrix collects the m eigenvectors spanning the Galactic signal subspace.

6.
We compute the effective dimension, m, of the foreground signal subspace (number of Galactic degrees of freedom, or principal components) by minimizing the Akaike Information Criterion (AIC, Akaike 1974): (5)where μ_{i} are the eigenvalues of .

7.
We apply the m–dimensional ILC filter (Remazeilles et al. 2011b) to the data in order to reconstruct the total Galactic signal at scale j: (6)where the estimated mixing matrix is given by (7)with collecting the m eigenvectors selected by the AIC criterion at scale j.
The reconstructed Galactic signal maps are finally synthesized as follows. We transform the estimated maps, , to spherical harmonic coefficients, then bandpass filter the harmonic coefficients by the respective needlet window, , and transform back to maps in real space. This operation provides one reconstructed Galactic signal map per needlet scale. We coadd these maps to obtain, for each frequency channel, the complete reconstructed Galactic signal map on the whole range of angular scales. The needlet windows are chosen so that , therefore conserving the total power in the synthesis.
The reconstruction of the CIB maps is performed as follows. In step 4, we replace Eq. (3) by so that the reconstructed signal is the sum of the Galactic dust plus the CIB. We then subtract the reconstructed Galactic dust (only), , from the Galactic dust plus CIB reconstruction.
It should be noted that the priors on the CMB and CIB angular power spectra are only used for estimating the dimension, m, of the Galactic signal subspace (step 5), not for the ILC filtering (step 7) in the reconstruction of the components of the emission.
We have validated the GNILC method on the Planck full focal plane simulations (Planck Collaboration XII 2016) before applying it to the Planck data. The results on simulations are presented in Sect. A.6 of Appendix A.
4. GNILC results on the thermal dust
Fig. 1 gnomonic projection of the sky centred at high latitude (l,b) = (90°,−80°). Top left: Planck 353GHz map. Top right: CMBremoved Planck 353GHz map. Middle left: dust model P13 at 353 GHz (MBB fit on CMBremoved Planck maps). Middle right: GNILC dust map at 353 GHz. Bottom left: dust model P15 at 545 GHz (Commander Bayesian fitting). Bottom right: GNILC dust map at 545 GHz. Maps at 353 GHz are shown at 5′ resolution, while maps at 545 GHz are smoothed to 7.5′ resolution. The GNILC dust maps have a nonuniform resolution (see Fig. 2) with 5′ resolution kept in regions of bright dust emission. In each image the local mean intensity has been subtracted for this comparison. 
Fig. 2 Effective beam FWHM of the GNILC dust maps on the whole sky (top panel) and on a area of the sky centred at high latitude (l,b) = (90°,−80°) (bottom panel). The spatially varying beam FWHM is the same for all frequencies. GNILC preserves the 5′scale power of the thermal dust in the high signaltonoise regions of the sky. 
Fig. 3 Fullsky map of the Galactic thermal dust emission: Planck 2013 (P13) thermal dust model at 353 GHz and 5′ resolution (top panel), suffering from CIB contamination at high latitudes, and the GNILC dust map (this work) at 353 GHz and 5′ resolution (bottom panel), for which the CIB is clearly filtered out at highlatitudes. A logarithmic colour scale is used here to highlight the lowintensity emission at high latitudes. The effective local beam of the GNILC dust maps is shown in Fig. 2. 
Fig. 4 GNILC dust maps at 353 GHz (left panel), 545 GHz (middle panel), and 857 GHz (right panel) on a gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). GNILC filters out the CIB anisotropies while preserving the smallscale dust signal (see bottom panel of Fig. 2). In these images, the local mean intensity of each map has been subtracted. 
Fig. 5 gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). Left: difference map (Planck 353 GHz −Planck CMB −GNILC dust) reveals the CIB anisotropies at 353 GHz. Middle: difference map (Planck 353 GHz – Planck CMB – dust model P13) revealing only the instrumental noise because the dust model P13, like the Planck observations at 353 GHz, still contains the CIB signal. Right: difference (dust model P13 – GNILC dust) revealing the amount of CIB contamination in the dust model P13 with respect to the GNILC dust map. In these images, the local mean intensity of each map has been subtracted. 
Fig. 6 Angular power spectra of the various maps at 353 GHz (top panel), 545 GHz (middle panel), and 857 GHz (bottom panel), on a fraction of the sky, f_{sky} = 57%: Planck map (solid black line), dust model P13 (dotted blue line, Planck Collaboration XI 2014), GNILC dust map (long dashed red line), GNILC dust map corrected for the residual noise (solid green line), and GNILC residual map (Planck map −Planck CMB −GNILC dust, solid yellow line), which is compared to the Planck CIB bestfit power spectrum (dashdot purple line, Planck Collaboration XXX 2014) and the Planck noise power spectrum (dashed orange line). 
4.1. Dust maps and power spectra
In Fig. 1 we compare various maps projected onto a high Galactic latitude area centred at (l,b) = (90°,−80°). In the top left panel of Fig. 1, the Planck353GHz channel map is shown. At 353 GHz the CMB radiation is clearly visible in the Planck observation map at high Galactic latitude, through degreescale temperature fluctuations which are typical in size of the CMB anisotropies. In the top right panel of Fig. 1, the Planck353GHz map is shown after subtraction of the Planck CMB map (i.e. the SMICA map from Planck Collaboration XII 2014). The CMBremoved Planck 353GHz map reveals the thermal dust emission, but is still quite noisy and contaminated by the CIB temperature anisotropies. The dust model P13 at 353 GHz, which has been computed by fitting an MBB spectrum to the CMBremoved Planck maps (Planck Collaboration XI 2014), is plotted in the middle left panel of Fig. 1. Because of the similar spectral signatures of the thermal dust and the CIB, the dust model P13 resulting from the spectral fitting can not avoid the leakage of the CIB fluctuations into the dust map. Conversely, in the GNILC dust map at 353 GHz produced in this work, the CIB anisotropies have been successfully filtered out, while the 5′scale dust emission has been conserved in the map. All the maps are shown at 5′ resolution, but the GNILC dust map has a local effective beam resolution that is shown in Fig. 2. The local beam resolution of the GNILC dust maps is not the result of a local smoothing of the maps, but the result of the thresholding of the needlet coefficients that depends on local signaltonoise ratio. In some highlatitude regions of the sky, beyond a certain angular scale, the power of the dust is found to be consistent with zero, i.e. the dimension of the Galactic signal subspace selected by the AIC criterion is m = 0 (see Fig. A.1) because the sky observations in this needlet domain become compatible with the CIBplusCMBplusnoise model.
In the bottom panels of Fig. 1, we compare the dust model P15 at 545 GHz from Planck Collaboration X (2016) and the GNILC dust map at 545 GHz produced in this work. The dust model P15 has been obtained by using a Bayesian fitting method, Commander (Eriksen et al. 2008), instead of the χ^{2} fitting method used for the dust model P13 in Planck Collaboration XI (2014). Commander makes an overall fit of many foreground parameters, including those of the thermal dust component (intensity, spectral index, and temperature). However, the Commander fitting again does not make any distinction between the CIB and the Galactic thermal dust, both sharing a similar MBB spectrum. In summary, thermal dust and CIB are still fitted as a single component in constructing the dust models P13 and P15. The dust model P15 at 545 GHz still shows CIB anisotropies at high latitude, whereas those CIB anisotropies have been successfully filtered out in the GNILC dust map at 545 GHz (see Fig. 1). Unlike the dust models P13 and P15, the GNILC dust maps are not the result of any fit of a dust model, but the result of a componentseparation procedure solely based on prior assumptions on the CIB, CMB, and noise angular power spectra.
Figure 3 shows the GNILC allsky map of the thermal dust at 353 GHz in the bottom panel. This map can be compared to the dust model P13 at 353 GHz from Planck Collaboration XI (2014), shown in the top panel. While the dust model P13 still shows visible smallscale contamination by CIB anisotropies at high latitude, in the GNILC dust map the CIB fluctuations are clearly filtered out at high latitude. The gnomonic projections, centred at (l,b) = (90°,−80°), of the various GNILC maps of the dust at 353, 545, and 857 GHz, are shown in Fig. 4.
At high Galactic latitude and small angular scales, the AIC criterion can select a dimension zero for the Galactic signal subspace, considering that in this region the Galactic signal is completely buried under the CIB and noise signals, and therefore the dust is compatible with zero. This aspect of the GNILC filtering is visible in the bottom panel of Fig. A.1, where in the highlatitude region at 5′ scale there are no Galactic degrees of freedom selected by the AIC criterion. Therefore, in practice the GNILC filtering is equivalent to a local smoothing of the sky map, depending on the relative power of the Galactic dust with respect to the local contamination by the CIB, the CMB, and the instrumental noise. The effective local beam FWHM over the sky of the GNILC dust maps is plotted in Fig. 2. Over 65% of the sky, where the dust signal is significant, the 5′ beam resolution is preserved by the GNILC filtering.
It is interesting to look at the residual map given by the difference between the CMBremoved Planck map and the dust map, i.e. the difference map (Planck map −Planck CMB map − dust map). In the case where the GNILC dust map is used in the subtraction, the residual map clearly shows the CIB anisotropies plus the instrumental noise (left panel of Fig. 5), as expected. Conversely, if the Planck 2103 dust model is used for the subtraction in place of the GNILC dust map then the residual map shows the instrumental noise only (middle panel of Fig. 5). This, again, indicates that the CIB anisotropies have leaked into the dust model P13, and therefore that the CIB anisotropies can not be recovered in the residual map. In the right panel of Fig. 5, we plot the difference between the Planck 2103 dust model and the GNILC dust map at 353 GHz, highlighting the amount of CIB leaking into the dust model P13.
The resulting angular power spectrum of the various dust maps and the residual maps at 353, 545, and 857 GHz are plotted in Fig. 6. We have used the HEALPix routine anafast (Górski et al. 2005) for computing the angular power spectrum of the maps. The amplitude of the power spectrum of the GNILC dust map at 353 GHz (long dashed red line) is reduced by a factor of 2 at ℓ ≈ 1000 with respect to the dust model P13 (dotted blue line) because of the removal of the CIB contamination. The power spectrum of the GNILC dust map is steeper than the power spectrum of the dust model P13. The GNILC dust power spectrum scales as a powerlaw C_{ℓ} ≈ ℓ^{2.7}. For completeness, the GNILC dust power spectrum corrected for residual noise is overplotted as a solid green line.
The angular power spectrum of the GNILC residual map, i.e. of the difference map (Planck map −Planck CMB −GNILC dust), typically shows the power of the CIB plus noise that has been filtered out in the GNILC dust map (solid yellow line). As illustrated in the top panel of Fig. 6, the angular power spectrum of the GNILC residual map successfully matches the sum of the Planck CIB bestfit power spectrum at 353 GHz (dashdot purple line) and the Planck 353GHz instrument noise power spectrum (dashed orange line). At 545 and 857 GHz the angular power spectrum of the GNILC residual map is below the Planck CIB bestfit power spectrum, which means that some amount of residual CIB contamination is still left in the GNILC dust maps at 545 and 857 GHz.
4.2. Thermal dust temperature and spectral index
Fig. 7 Fullsky thermal dust parameter maps: temperature (top row), spectral index (middle row), and map of the χ^{2} statistic of the fit (bottom row). Left panels: GNILC modified blackbody (MBB) fit. Right panels: PR2 modified blackbody (MBB) fit a la model P13. 
Following Planck Collaboration XI (2014), we fit in each pixel a modified blackbody (MBB) spectral model to the GNILC dust maps at 353, 545, 857, and 3000 GHz in order to estimate the dust temperature, spectral index, and optical depth over the sky. We also performed an analysis similar to that of Planck Collaboration XI (2014), which used the PR1 data, by fitting the same MBB model to the CMBremoved PR2 maps, in place of the GNILC dust maps. We refer to this as the PR2 MBB fit. This allows us to highlight the improvement in estimating the dust temperature and spectral index after filtering out the CIB anisotropies with GNILC.
Fig. 8 gnomonic projections of the dust temperature maps at high latitude, b = −80° (top panels) and low latitude, b = −20° (bottom panels). Left: GNILC MBB fit. Right: PR2 MBB fit according to model P13. 
Fig. 9 gnomonic projections of the dust spectral index maps at high latitude, b = −80° (top panels) and low latitude, b = −20° (bottom panels). Left: GNILC MBB fit. Right: PR2 MBB fit according to model P13. 
Fig. 10 Highlatitude area of the sky with f_{sky} = 21% (top) and lowlatitude area of the sky with f_{sky} = 20% (bottom) that are considered in Fig. 11 and Table 1. 
4.2.1. χ^{2} fitting
The model of dust emission that we fit to the data is a modified blackbody (MBB) spectrum with three parameters: (8)where ν_{0} = 353 GHz is the reference frequency, τ_{0}(p) the dust optical depth at 353 GHz in pixel p, T(p) the dust temperature in pixel p, and β(p) the dust spectral index in pixel p. The function is the Planck law for blackbody radiation.
We use a standard χ^{2} fitting method as in Planck Collaboration XI (2014). However, there a twostep approach was adopted for the fit; in order to limit the fluctuations in the estimated parameters induced by the noise and the CIB contamination, the spectral index parameter was estimated at 30′ resolution in a first step, then the temperature and the optical depth were fit at 5′ resolution. Given that we already have cleaned the thermal dust from CIB contamination at 353, 545, 857, and 3000 GHz by using the GNILC componentseparation method, there is no reason to perform a lowresolution MBB fit on the cleaned GNILC maps. Therefore, we will fit the three parameters τ_{0}, β, and T simultaneously at full resolution (5′), instead of following the twostep approach adopted in Planck Collaboration XI (2014). For the fit, we use the frequency data at 353, 545, 857, and 3000 GHz, either from the unfiltered PR2 data (i.e. inputs similar to those used to produce the dust model P13) or from the CIBremoved GNILC dust maps. In this way, we will highlight the improvement in the estimated dust parameters resulting from the removal of the CIB fluctuations with GNILC.
In most of the images presented in this paper the local average of the dust maps has been subtracted to facilitate sidebyside comparisons of the different versions of the Planck dust map in terms of contamination by CIB fluctuations. There is no subtraction of any offset in the released GNILC products themselves. In order to fit for the dust spectral parameters, τ_{0}, T, and β, the offsets of the GNILC dust maps have been estimated by correlation with the H i map at high latitude in the exact same way as described in Planck Collaboration XI (2014). The offsets of the GNILC dust maps are found to be 0.1248, 0.3356, 0.5561, and 0.1128 MJy sr^{1} at 353, 545, 857, and 3000 GHz respectively. The uncertainties on the absolute calibration of the Planck channels have been estimated from the observation of planets (Planck Collaboration VIII 2016); they are 1.2% at 353 GHz, 6.3% at 545 GHz, and 6.1% at 857 GHz. The calibration incertainty at 3000 GHz is 13.5% (MivilleDeschênes & Lagache 2005). Calibration uncertainties and offset uncertainties have been taken into account in the χ^{2} fitting, following the procedure detailed in the appendix of Planck Collaboration XI (2014).
4.2.2. Parameter maps
The results of the MBB fit to the GNILC dust maps are shown on the left panels of Fig. 7. The estimated GNILC temperature map and GNILC spectral index map are compared to the PR2 MBB fit temperature map and the PR2 MBB fit spectral index map at 5′ resolution. The PR2 MBB fit is similar to the dust model P13 of Planck Collaboration XI (2014), i.e. the CIB anisotropies have not been filtered out, except that the model fitting is applied to the PR2 data instead of the PR1 data and not performed in two steps, but carried out simultaneously for the three dust parameters at 5′ resolution. The impact of the CIB contamination on the measurement of the dust temperature and dust spectral index is particularly significant at high latitude in the PR2 MBB fit.
In the bottom panels of Fig. 7, we plot the resulting χ^{2} map of both the GNILC MBB fit and the PR2 MBB fit. This provides a direct measurement of the goodnessoffit. The reason for some reduced χ^{2} values being smaller than unity is mostly that calibration uncertainties are included per pixel in the fit, to give less weight to data points with larger uncertainty (Planck Collaboration XI 2014). However, the exact scale of χ^{2} is not relevant here, what is important is that the pixeltopixel differences in the goodnessoffit are strongly reduced for the GNILC MBB fit because of the removal of the CIB contamination in the GNILC dust maps. Clearly, the CIBfiltered GNILC maps lead to a better fitting of the MBB model over the sky than the unfiltered PR2 maps. Near the Galactic plane, the χ^{2} values between the GNILC MBB fit and the unfiltered PR2 MBB fit are consistent because the CIB contamination plays a negligible role where the dust emission is bright. For a given spectral model of thermal dust, here a single MBB model, the GNILC maps provide higher precision than the unfiltered PR2 maps because of the removal of CIB fluctuations prior to fitting. However, the MBB model might not be the best parametrization of the thermal dust emission in the inner Galactic plane region, which shows high values of the χ^{2} statistic for both PR2 and GNILC fits. It is likely that multiple MBB components of dust might contribute to the emission along the line of sight, in which case the exact parametrization of the thermal dust spectral energy distribution in the inner Galactic plane region is not trivial and might need more than three effective parameters.
Mean and dispersion of the dust temperature and spectral index in different areas of the sky (full sky, high latitude, low latitude). Top: GNILC MBB fit. Middle: PR2 MBB fit a la dust model P13. Bottom: dust model P15 (Commander60′).
Fig. 11 Normalized histograms of T_{dust} and β_{dust} at 5′ resolution for the GNILC MBB fit (red contours) and the PR2 MBB fit a la model P13 (green contours). The normalized histograms for the dust model P15 (Commander fit at 60′ resolution) are overplotted (blue contours). The histograms are computed from the subset of pixels corresponding to either the highlatitude area in the sky with f_{sky} = 21% (upper panels), the lowlatitude area in the sky with f_{sky} = 20% (middle panels), or the whole sky (lower panels). Due to CIB contamination at highlatitude, the PR2 MBB fits show larger dispersion than the GNILC MBB fits in the distributions of T_{dust} and β_{dust}. 
Fig. 12 GNILC CIB maps for a large fraction of the sky at 353 GHz (top), 545 GHz (middle), and 857 GHz (bottom). Apart from thermal dust reconstruction, the GNILC componentseparation method gives access to CIB anisotropies over 57% of the sky. The left panels show the GNILC CIB maps at full resolution while the right panels show the CIB maps smoothed to one degree resolution. 
Fig. 13 GNILC CIB maps at 353 GHz (left), 545 GHz (middle), and 857 GHz (right) on a gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). The partial spatial correlation of the CIB anisotropies between pairs of frequencies is clearly visible. In this figure, the local mean intensity of each map has been subtracted. 
In Figs. 8 and 9 we compare the GNILC and PR2 MBB fit for temperature and spectral index respectively, at low and high latitudes in the sky. The improvement from GNILC in terms of the reduction of the CIB contamination is particularly visible at high latitude for both temperature and spectral index.
In Fig. 10 we define low and highlatitude areas of the sky to look at the evolution of the distribution of the dust temperature and spectral index with respect to latitude. Figure 11 shows the normalized histograms of the temperature map, T, and the spectral index map, β, for three different fits: GNILC (red contours); PR2 MBB fit a la model P13 (green contours); and dust model P15 (blue contours). It is important to note that the dust model P15 is a lowresolution Bayesian fit at 60′ resolution with Gaussian priors on T (23 ± 3 K) and β (1.55 ± 0.1). We distinguish three areas in the sky: the highlatitude area defined in Fig. 10, covering 21% of the sky (top panels); the lowlatitude area defined in Fig. 10, covering 20% of the sky (middle panels); and the whole sky (bottom panels). As a complement to Fig. 11, Table 1 summarizes the mean bestfit values of the dust parameters, along with their 1σ errors, for the three different products (GNILC MBB fit, PR2 MBB fit similar to the dust model P13, and dust model P15) in the three different areas of the sky. The histograms in Fig. 11 highlight the impact of the CIB anisotropies on the dust spectral parameters: at high latitude the CIB contamination increases the scatter in the temperature and spectral index distributions. The removal of the CIB anisotropies with GNILC reduces the dispersion in the dust temperature by 40% at high latitude (30% on the whole sky) with respect to the PR2 MBB fit and by 10% with respect to the dust model P15 (Table 1), even though the P15 temperature fit is smoothed to 60′ resolution. The impact of the CIB removal with GNILC is even more significant for the dust spectral index, with the dispersion reduced by 60% at high latitude (50% on the whole sky) with respect to the PR2 MBB fit. The 1σ error on the P15 spectral index has a lowest value of 0.05 in any area of the sky for two reasons: first, the resolution of the P15 spectral index is much lower (60′); second, a tight prior has been imposed on β in the Commander fit (Planck Collaboration X 2016). With GNILC we find a dust temperature of T = (19.4 ± 1.3) K and a dust spectral index of β = 1.6 ± 0.1 as the bestfit values on the whole sky, where the error bars show the dispersion over the sky of the parameter values.
5. GNILC results on the CIB
5.1. CIB maps
Fig. 14 Maps of cosmic infrared background (CIB) anisotropies in GHIGLS fields (Martin et al. 2015). Left column: CIB 2013 (Planck Collaboration XXX 2014). Middle column: CIB GNILC. Right column: difference (CIB 2013 − GNILC). From top to bottom row, the GHIGLS fields are: Boötes; MC; N1; and SP. The size and the location in the sky of each field are defined in Table 1 of Martin et al. (2015). 
Fig. 15 T–T scatter plot between CIB GNILC and CIB 2013. The fields are Boötes, MC, N1, and SP. 
The GNILC method is flexible by allowing either the recovery of the dust map with the CIBplusCMBplusnoise filtered out (shown in this paper) or the recovery of the dustplusCIB map with the removal of the CMBplusnoise only, depending on whether or not one uses the prior on the CIB power spectrum. Therefore, from the difference between the unfiltered GNILC dustplusCIB map and the CIBfiltered GNILC dust map we are able to reveal the CIB anisotropies at different frequencies over a large area of the sky. The resulting GNILC CIB maps at 353, 545, and 857 GHz, reconstructed over a large fraction of the sky, are shown in Fig. 12, both at full resolution and smoothed to 1° resolution. We can see residual zodiacal light emission along the ecliptic plane in the lowresolution CIB maps. This residual comes from the combination by GNILC of different data sets, namely the Planck, IRAS, and SFD maps, in which the zodiacal light has been corrected differently. The zodiacal light emission has been reduced in Planck data to a negligible level compared to the CMB and dust emissions but is still at a level comparable to the amplitude of the CIB emission.
The GNILC products give us access to the CIB anisotropies on a much larger fraction of the sky (approximately 57%) than the Planck CIB maps produced in Planck Collaboration XXX (2014). The GNILC CIB maps can be used as tracers of the dark matter distribution because they allow for an exploration over large areas of the sky of the crosscorrelations between CIB and CMB lensing fields (Planck Collaboration XVIII 2014) or between CIB and other tracers of largescale structure (Serra et al. 2014).
In Fig. 13 we show the GNILC CIB maps at 353, 545, and 857 GHz in a highlatitude region of the sky centred at (l,b) = (90°,−80°). The partial spatial correlation of the CIB anisotropies between pairs of frequencies decreases when the ratio between the two frequencies is further from unity, as expected from the redshift distribution of the CIB anisotropies (Planck Collaboration XXX 2014).
5.2. GNILC CIB versus CIB 2013 in small fields
We now check the consistency between the new GNILC CIB maps and the CIB 2013 maps from Planck Collaboration XXX (2014). In Fig. 14, we compare the GNILC CIB map at 545 GHz with the CIB 2013 map at 545 GHz in four different fields of the Green Bank Telescope H i Intermediate Galactic Latitude Survey (GHIGLS) defined in Martin et al. (2015): Boötes, MC, N1, and SP. The difference map (CIB 2013 −GNILC CIB) is also shown within the same fields. In Fig. 15 we plot the T–T correlation between the GNILC CIB map and the CIB 2013 map in the common fields. We have used a leastsquares bisector linear regression (Isobe et al. 1990) for computing the correlation coefficient between both products.
Figures 14 and 15 show that the GNILC CIB maps are consistent with the CIB 2013 maps within the fields considered. In particular the Pearson correlation coefficient between the GNILC CIB maps and the CIB 2013 maps is larger than 0.8 in all fields, with a T–T slope of 0.998 ± 0.005 in the Boötes field, 0.958 ± 0.006 in the MC field, 0.972 ± 0.006 in the N1 field, and 0.935 ± 0.006 in the SP field. Despite the high correlation between the GNILC CIB maps and the CIB 2013 maps within the GHIGLS fields, the correlation is not perfect because both sets of maps were produced from different data releases, respectively PR2 and PR1, for which there have been changes in the calibration coefficients. In addition, estimates of the optical beam resolution have slightly changed between both data releases, therefore not guaranteeing the exact same resolution of both CIB products.
6. Correlations of the CIB and dust maps with the H I map
In the diffuse interstellar medium the dust emission is tightly correlated with the line emission of neutral hydrogen (H i). In this respect, as a tracer of the Galactic dust emission the H i emission map can be used to detect any residual Galactic dust emission in the GNILC CIB maps.
We compute the T–T correlation between the H i map (local plus intermediate velocity clouds) of the LAB survey (Kalberla et al. 2005; Land & Slosar 2007) and the GNILC CIB maps in order to detect any Galactic residual in the reconstructed CIB maps. Figure 16 shows that the correlation between the GNILC CIB map at 353 GHz and the H i map is consistent with zero (Pearson correlation coefficient of 0.004), therefore showing no significant residual Galactic emission in the GNILC CIB map.
Fig. 16 T–T scatter plot between the H i map and the GNILC CIB map at 353 GHz over 57% of the sky (see Fig. 12). 
We also compute the T–T correlation between the dust map at 353 GHz and the H i map at high Galactic latitude. The high latitude region is defined as the area of the sky where the local beam FWHM of the GNILC dust map is larger than 15′ (Fig. 2). The maps are degraded to HEALPix N_{side} = 256.
Fig. 17 T–T scatter plot at high Galactic latitude between the H i map and the dust optical depth maps at 353 GHz: dust model P13 (green), PR2 MBB fit (black), GNILC MBB fit (blue). While the slope of the correlation with the H i map is consistent for all the dust optical depth maps, the scatter is smallest for the GNILC dust optical depth map because of the removal of the CIB temperature anisotropies. 
Figure 17 shows the scatter plot between the H i map and the dust optical depth map for the dust model P13 (green), the PR2 MBB fit (black), and the CIBfiltered GNILC MBB fit (blue). The PR2 MBB fit shows larger scatter than the dust model P13 because in the former the fit is performed in one step at 5′ resolution for all three dust parameters, while in the latter the fit was performed in two steps, with the spectral index first fitted at lower resolution (30′), which slightly reduces the scatter due to CIB contamination. Owing to the large reduction of the CIB contamination in the GNILC maps, the correlation with the H i shows the most reduced scatter.
Fig. 18 Polar orthographic projections of the E(B−V) maps at north (left) and south (right) Galactic poles in the PanSTARRS 1 survey area (Onaka et al. 2008). Top: E(B−V) map from Green et al. (2015). Bottom: GNILC E(B−V) map. 
Fig. 19 T–T scatter plot between the ratio E(B−V) /N_{H} and the gas column density N_{H} for the twodimensional projection of the E(B−V) map of Green et al. (2015) (red diamonds) and the GNILC E(B−V) map (blue diamonds) in the highlatitude region of the sky defined in Fig. 10. Each point is the average of E(B−V) /N_{H} values in a bin of N_{H}. The bin size varies such that there is always the same number of samples per bin. 
The particles producing thermal dust emission also cause extinction of the light from stars and quasars, which is quantified by Galactic reddening, E(B−V). We develop a GNILC reddening E(B−V) map by multiplying the GNILC dust optical depth map, τ_{353}, by the factor 1.49 × 10^{4} mag derived in Planck Collaboration XI (2014) from the correlation between the reddening of quasars and dust optical depth along the same line of sight. Recently, Green et al. (2015) have produced a threedimensional dust reddening E(B−V) map based on stars in the PanSTARRS 1 survey. We projected their threedimensional dust reddening into a twodimensional map by computing the median reddening in the farthest distance bin. Figure 18 compares this with the GNILC E(B−V) map, focusing on common area in the north and south Galactic caps where the gas column density is low. Although it shows spatial structure similar to the GNILC E(B−V) map, the E(B−V) map from Green et al. (2015) is noisier. Just as there is a good correlation of τ_{353} and the H i gas column density, N_{H}, the correlation of reddening E(B−V) and N_{H} has been long established (e.g., Savage & Jenkins 1972). To explore this further for the low column density area of the sky defined in Fig. 10, in Fig. 19 we plot the ratio E(B−V) /N_{H} binned with respect to N_{H} for both GNILC and the reddening map from Green et al. (2015). For GNILC the trend with N_{H} is quite flat; the horizontal line plotted is compatible with the behaviour for the opacity τ_{353}/N_{H} found in Planck Collaboration XI (2014) over the same range in N_{H} (see their Fig. 20). On the other hand, for the E(B−V) map from Green et al. (2015) we see a strong dependence of the ratio on N_{H}, i.e. a lack of linearity between this measure of E(B−V) and N_{H}. The binned results from the map of Green et al. (2015) also show a larger dispersion. From these three perspectives, the GNILC optical depth map appears to provide a better template for E(B−V) studies. Both products show an excess ratio at the lowest column densities (N_{H}< 1 × 10^{20} cm^{2}); this could be the signature of dust mixed with ionized hydrogen, which is not traced by neutral H i emission (see discussion in Planck Collaboration XI 2014).
7. Conclusions
We have produced a significantly improved allsky map of the Galactic thermal dust emission from the Planck data and the IRAS 100 μm map. We have fitted a modified blackbody model in each pixel to the GNILC dust maps produced at 353, 545, 857, and 3000 GHz. The new PlanckGNILC dust model has been compared with the dust models P13 and P15 and has been shown to be significantly less contaminated by CIB and noise.
By exploiting the distinct signature of Galactic dust and extragalactic CIB angular power spectra, the GNILC method successfully separates the Galactic thermal dust emission from the CIB anisotropies in the Planck PR2 maps. We have reduced the dispersion due to CIB contamination of the estimated dust temperature and spectral index in the GNILC dust map with respect to the P13 dust map by a factor 1.3 for T on the whole sky (1.6 at high latitude) and a factor 2.1 for β on the whole sky (2.6 at high latitude). The GNILC dust map at 353 GHz has already been implemented as the thermal dust model in the released Planck simulations of the sky (Planck Collaboration XII 2016).
The GNILC method, presented in this work, also gives access to the CIB anisotropies over a large fraction of the sky. Within the small fields considered in Planck Collaboration XXX (2014), the GNILC CIB maps and the CIB 2013 maps are found to be consistent with a Pearson correlation coefficient larger than 0.8. The GNILC CIB maps can be very useful as indirect tracers of the dark matter over a large area of the sky and they are recommended for the investigation of crosscorrelations with galaxy weak lensing data and other tracers of largescale structure.
The new PlanckGNILC products are made publicly available on the Planck Legacy Archive^{3}. They include:

the CIBremoved GNILC thermal dust maps at 353, 545, and 857 GHz;

the GNILC CIB maps at 353, 545, and 857 GHz;

the GNILC dust optical depth map;

the GNILC dust spectral index map;

the GNILC dust temperature map;

the GNILC effective beam map for the thermal dust.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
http://pla.esac.esa.int/pla; see Release PR2 – 2015, Foreground maps, CIB or Dust.
It could be argued that the nonGaussianity of the CIB anisotropies should be taken into account in the statistics, in order to make the separation of Galactic foregrounds and CIB more accurate. However, the nonGaussianity of the CIB is negligible compared to the nonGaussianity of the Galactic foregrounds, so that Gaussian statistics is a sufficient approximation for separating the thermal dust and the CIB.
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planckcollaboration. Some of the results in this paper have been derived using the HEALPix package. The research leading to these results has received funding from the ERC Grant No. 307209.
References
 Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716 [Google Scholar]
 Basak, S., & Delabrouille, J. 2012, MNRAS, 419, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Basak, S., & Delabrouille, J. 2013, MNRAS, 435, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Hill, R. S., Hinshaw, G., et al. 2003, ApJS, 148, 97 [NASA ADS] [CrossRef] [Google Scholar]
 BICEP2/Keck Array and Planck Collaborations 2015, Phys. Rev. Lett., 114, 101301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Cardoso, J., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 735 [Google Scholar]
 Delabrouille, J., Cardoso, J.F., & Patanchon, G. 2003, MNRAS, 346, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Draine, B. T., & Li, A. 2007, ApJ, 657, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1 [NASA ADS] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Guilloux, F., Faÿ, G., & Cardoso, J.F. 2009, Appl. Comput. Harmon. Anal., 26, 143 [CrossRef] [Google Scholar]
 Isobe, T., Feigelson, E. D., Akritas, M. G., & Babu, G. J. 1990, ApJ, 364, 104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalberla, P. M. W., Burton, W. B., Hartmann, D., et al. 2005, A&A, 440, 775 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kelsall, T., Weiland, J. L., Franz, B. A., et al. 1998, ApJ, 508, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Kullback, S. 1968, Information theory and statistics (New York: Dover) [Google Scholar]
 Lagache, G., Puget, J.L., & Dole, H. 2005, ARA&A, 43, 727 [NASA ADS] [CrossRef] [Google Scholar]
 Land, K., & Slosar, A. 2007, Phys. Rev. D, 76, 087301 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153 [NASA ADS] [CrossRef] [Google Scholar]
 MivilleDeschênes, M.A., & Lagache, G. 2005, ApJS, 157, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Narcowich, F., Petrushev, P., & Ward, J. 2006, SIAM J. Math. Anal., 38, 574 [CrossRef] [Google Scholar]
 Olivari, L. C., Remazeilles, M., & Dickinson, C. 2016, MNRAS, 456, 2749 [NASA ADS] [CrossRef] [Google Scholar]
 Onaka, P., Tonry, J. L., Isani, S., et al. 2008, Proc. SPIE, 7014, 70140D [CrossRef] [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2011, A&A, 536, A19 [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 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 XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2016, A&A, 594, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2016, A&A, 594, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2016, A&A, 594, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2016, A&A, 594, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2016, A&A, 594, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XVII. 2014, A&A, 566, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXVIII. 2015, A&A, 582, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXX. 2016, A&A, 586, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XXXII. 2016, A&A, 586, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puget, J.L., Abergel, A., Bernard, J.P., et al. 1996, A&A, 308, L5 [NASA ADS] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011a, MNRAS, 410, 2481 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Delabrouille, J., & Cardoso, J.F. 2011b, MNRAS, 418, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Aghanim, N., & Douspis, M. 2013, MNRAS, 430, 370 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Dickinson, C., Banday, A. J., BigotSazy, M.A., & Ghosh, T. 2015, MNRAS, 451, 4311 [NASA ADS] [CrossRef] [Google Scholar]
 Remazeilles, M., Dickinson, C., Eriksen, H. K. K., & Wehus, I. K. 2016, MNRAS, 458, 2032 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Press, W. H., & Kirshner, R. P. 1996, ApJ, 473, 588 [NASA ADS] [CrossRef] [Google Scholar]
 Savage, B. D., & Jenkins, E. B. 1972, ApJ, 172, 491 [NASA ADS] [CrossRef] [Google Scholar]
 Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Serra, P., Lagache, G., Doré, O., Pullen, A., & White, M. 2014, A&A, 570, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Description of the GNILC method
The generalized needlet internal linear combination, GNILC (Remazeilles et al. 2011b), is a componentseparation method designed to reconstruct the diffuse emission of a complex component originating from multiple correlated sources of emission, such as the Galactic foreground emission or the cosmic infrared background radiation.
GNILC is a multidimensional generalization (Sect. A.1) of the standard internal linear combination (ILC) method, which has been extensively used to reconstruct onedimensional components such as the CMB emission (Bennett et al. 2003; Planck Collaboration XII 2014) or the SunyaevZeldovich (SZ) signal (Remazeilles et al. 2013; Planck Collaboration XXII 2016). A twodimensional extension of the ILC, the socalled Constrained ILC, was first developed by Remazeilles et al. (2011a) to reconstruct both maps of the CMB and the SZ components, with vanishing contamination from one into the other. GNILC is a further generalization in which the dimension of the ILC filter, which is related to the dimension of the signal subspace, is no longer fixed, but varies with both the direction in the sky and the angular scale, depending on the effective local signaltonoise ratio, i.e. the local conditions of contamination in both real space and harmonic space.
In this work, the signal is Galactic and the noise contributions consist of the CIB, the CMB, and the instrumental noise. The effective signaltonoise ratio is determined locally both over the sky and over different angular scales by decomposing the data onto a wavelet (needlet) frame and by making use of a prior on the CIB power spectrum (Sect. A.2). The dimension of the Galactic signal is estimated locally in each wavelet domain through a modified principal component analysis (PCA) constrained by the power spectrum of the CIB (Sect. A.3). The effective dimension of the Galactic subspace, given by the local number of principal components, is not determined ad hoc but through a statistical model selection by the Akaike Information Criterion (Sect. A.4). The prior of the CIB power spectrum is only used at the stage of determining the number of principal components, i.e. the dimension of the Galactic signal subspace. There is no prior assumption about the Galactic signal.
Appendix A.1: Multidimensional internal linear combination
We model the sky observation, x_{i}(p), at frequency channel i and in the direction p in the sky (pixel), as the combination of the Galactic foreground emission, the CIB emission, the CMB emission, and the instrumental noise: (A.1)assuming no correlations between the different components of emission. Equation (A.1) can be recast in the N_{ch} × 1 vector form, where N_{ch} is the number of frequency channels: (A.2)Here x = (x_{i})_{1 ≤ i ≤ Nch} collects the N_{ch} observation maps, each of them being a mixture of the Galactic foreground emission, f (i.e. the thermal dust emission at high frequencies), the CIB emission, s^{CIB}(p), the CMB emission, s^{CMB}(p), scaling with a known spectral distribution, a, and the instrumental noise, n.
The N_{ch} × N_{ch} frequencyfrequency covariance matrix of the sky observations, R(p) = (R_{ij}(p))_{1 ≤ i,j ≤ Nch} = ⟨ x(p)x(p)^{T} ⟩, is (A.3)where R_{f} = ⟨ ff^{T} ⟩ is the covariance matrix of the Galactic signal, the covariance matrix of the CIB, the covariance matrix of the CMB, and R_{noise} = ⟨ nn^{T} ⟩ the covariance matrix of the noise.
The Galactic foreground signal, f, is a complex multicomponent emission emanating from various physical processes (e.g., thermal dust emission, synchrotron emission, and freefree emission) with spectral properties varying over the sky. The number of degrees of freedom of the Galactic foreground signal would be infinite in the case of an infinitely narrow beam. However, in practice the beam is finite and the observations are limited by the number of frequency channels and the level of noise so that the effective number of Galactic degrees of freedom required to describe the Galactic emission is finite. In addition, the various physical components of the diffuse Galactic emission are correlated, therefore the Galactic signal, f, can be represented as the superposition of a relatively small number, m, of independent (not physical) templates, t: (A.4)where F is an N_{ch} × m mixing matrix giving the contribution from the templates in each frequency channel. Therefore, the covariance matrix of the Galactic signal is an N_{ch} × N_{ch} matrix of rank m: (A.5)where R_{t} = ⟨ tt^{T} ⟩ is a fullrank m × m matrix.
We now address the problem of estimating the set of maps, f(p), i.e. determining a “catchall” foreground component comprising the emission of the diffuse Galactic interstellar medium. The objective is to construct estimated maps, , which are good matches to what would be observed by the instrument in the absence of CMB, CIB, and noise.
For extracting such an emission component from multifrequency observations, we propose to generalize the internal linear combination (ILC) method to address the case of such a “multidimensional component” (m–dimensional, with m<N_{ch}). We consider the estimation of f by a weighted linear operation (A.6)where the N_{ch} × N_{ch} weight matrix, W, is designed to offer unit response to the Galactic foreground emission while minimizing the total variance of the vector estimate . Stated otherwise, the matrix W is the minimizer of E(   Wx   ^{2}) under the constraint WF = F. The weights matrix, W, thus solves the following constrained variance minimization problem (A.7)where is the covariance matrix of the observations, x. This problem can be solved by introducing a Lagrange multiplier matrix, Λ, and the Lagrangian (A.8)By differentiating Eq. (A.8) with respect to W, one finds that ∂ℒ(W,Λ) /∂W = 0 is solved by (A.9)By imposing the constraint WF = F on Eq. (A.9), one then finds that . Hence, the solution of the Eq. (A.7) is given by the ILC weight matrix (A.10)Multidimensional ILC appears as a direct generalization of the onedimensional ILC of Bennett et al. (2003). The mixing matrix, F, of the Galactic signal and its dimension, m, are the unknowns of the problem. However, it is important to notice that expression (A.10) for W is invariant if F is changed into FT for any invertible matrix T. Hence, implementing the multidimensional ILC filter (A.10) only requires that the foreground mixing matrix, F, be known up to right multiplication by an invertible factor (Remazeilles et al. 2011b). In other words, the only meaningful and mandatory quantity for implementing a multidimensional ILC is not the “true” mixing matrix but the column space of F, i.e. the dimension, m, of the Galactic signal subspace. This is the purpose of the Sects. A.3 and A.4.
The frequencyfrequency covariance matrix of the observations, R, can be estimated from the observation maps across the frequency channels. The coefficients of the covariance matrix of the observations for the pair of frequencies (a,b) is computed empirically as follows: (A.11)where is a domain of pixels centred around the pixel p. In this work, the pixel domain, , is defined by smoothing the product map x_{i}(p) x_{j}(p) with a Gaussian window in pixel space.
Appendix A.2: Nuisance covariance matrix
We define the “nuisance” covariance matrix, R_{N}, as the sum of the CIB covariance matrix, the CMB covariance matrix, and the noise covariance matrix: (A.12)Therefore, (A.13)Our aim is to obtain an estimate, , of the nuisance covariance matrix. In combination with the estimate of the full covariance matrix, (Eq. (A.11)), the nuisance covariance matrix, , allows us to constrain the Galactic signal subspace.
The noise covariance matrix, R_{noise}, can be estimated from the covariance of the halfdifference maps between the HM1 and the HM2 halfmission surveys of Planck (Planck Collaboration VIII 2016), since the sky emission cancels out in the difference while the noise does not. The halfdifference map for frequency channel a is (A.14)where and are the Planck HM1 map and the Planck HM2 map, respectively. Similarly to the full covariance matrix in Eq. (A.11), the noise covariance matrix is then estimated empirically: (A.15)The Planck CMB bestfit ΛCDM model (Planck Collaboration XV 2014) and the Planck CIB bestfit models (Planck Collaboration XXX 2014) are used as priors to estimate: (i) the CMB covariance matrix, R_{CMB}; and (ii) the CIB covariance matrix, R_{CIB}, in pixel space.
We first simulate a Gaussian CMB map, , having a power spectrum given by the Planck CMB bestfit , and we scale it across frequencies through the known CMB spectral distribution, a. In harmonic space, the crosspower spectrum of the simulated CMB maps is given by (A.16)where is the Planck CMB bestfit. From the simulated CMB maps, , we are able to compute the CMB covariance matrix in pixel space, , in the same way as in Eqs. (A.11) and (A.15).
From the Planck CIB bestfit crossand autopower spectra, , we simulate N_{ch} correlated Gaussian maps^{4}, , having a correlation matrix in harmonic space given by (A.17)where are the Planck CIB bestfits. From the simulated correlated CIB maps, , we can compute the CIB covariance matrix in pixel space, , in the same way as in Eqs. (A.11) and (A.15). For all a,b< 143 GHz, the coefficients of are set to zero because there is no Planck measurement of the CIB at those low frequencies, and therefore they are considered negligible.
Appendix A.3: Determining the Galactic signal subspace with a constrained PCA
Once both the observation covariance matrix, , and the nuisance covariance matrix, , have been computed, we can “whiten” the Planck data by making the transformation (A.18)such that the covariance matrix of the transformed Planck observations is now given by (A.19)Assuming that the prior CIB covariance matrix, , is close the real CIB covariance matrix, R_{CIB}, we have that (A.20)where I is the identity matrix. In this way, from Eq. (A.13) the covariance matrix of the transformed Planck observations becomes (A.21)such that the power of all the nuisance contamination, including CIB, CMB, and noise, is encoded in the matrix I, which is close to an identity matrix. Typically, the coefficients of the transformed observation covariance matrix (Eq. (A.21)) provide the signaltonoise ratio over the sky, i.e. the power of the Galactic signal divided by the overall power of the CIBplusCMBplusnoise.
Therefore, by diagonalizing the transformed observation covariance matrix (Eq. (A.21)), we obtain the following eigenstructure:
(A.22)In such a representation, the eigenvalues of the covariance matrix that are close to unity do not contain any relevant power of the Galactic signal and the signal is dominated by the CIB, the CMB, and the noise. The corresponding eigenvectors span the nuisance subspace characterized by the number of degrees of freedom, N_{d.o.f.}, of the CIB plus CMB plus noise only. Conversely, the subset of eigenvectors collected in the N_{ch} × m matrix U_{S}, for which the eigenvalues of significantly depart from unity, span the Galactic signal subspace (principal components). The number of eigenvalues, m, that are much larger than unity corresponds to the dimension of the Galactic signal subspace, i.e. the number of independent (unphysical) templates contributing to the Galactic signal. This is a constrained principal component analysis (PCA), in the sense that the PCA is driven by the local signaltonoise ratio.
The diagonalization of the transformed Planck covariance matrix of Eq. (A.22) can be written in a compact form as (A.23)where (A.24)is a m × m diagonal matrix and [U_{S}  U_{N}] is an N_{ch} × N_{ch} orthonormal matrix collecting all the eigenvectors of . Using Eq. (A.23) and the orthonormality condition, , the covariance matrix of the Galactic signal can be written as (A.25)This is the power that “best matches” what would be observed by the instrument in the absence of CIB, CMB, and noise. Therefore, we estimate the Galactic signal, f, by (A.26)where (A.27)is an N_{ch} × m mixing matrix estimate, and t is an m × 1 vector of independent templates whose covariance matrix is given by (A.28)The estimate, , is the only useful information for implementing the multidimensional ILC filter of Eq. (A.10). It can be different from the true mixing matrix, F, of the Galactic signal as long as the column space is the same: let T be some invertible m × m matrix and consider the transformed matrices and . These transformed matrices are an alternate factorization of the covariance matrix of the Galactic signal, , but they are equivalent, since by construction, . The ILC weights of Eq. (A.10) are unchanged under right multiplication by an invertible matrix. Therefore, the m independent templates, t, can be replaced by any other linear combination, Tt, as long as we are interested in reconstructing the overall Galactic signal, f.
The number, m, of principal components in Eq. (A.22), or the effective dimension of the Galactic signal subspace, can vary over the sky, depending on the local signaltonoise ratio. In particular, at high Galactic latitudes this number decreases because the contributions of the CIB, the CMB, and the noise starts dominating the Galactic signal (Fig. A.1). The dimension of the Galactic signal subspace is also expected to vary with angular scale because at small angular scales the power of the Galactic signal becomes dominated by the power of the noise and the CIB. Therefore, for the accurate separation of the Galactic and CIB signals we find it useful to estimate the dimension of the Galactic subspace locally both in space and in scale. This is achieved by decomposing the data on a wavelet frame.
In this work, the analysis is performed on a needlet frame (see, e.g., Delabrouille et al. 2009; Remazeilles et al. 2013 for the use of needlets in component separation). Basically, the spherical harmonic transforms of the maps, a_{ℓm}, are bandpass filtered in harmonic space, then transformed back into real space, therefore conserving a specific range of angular scales in the map. The result is called a needlet map, characterized by a given range of angular scales. The multidimensional ILC Eq. (A.6) is performed on the needlet maps independently for each range of scale, and the synthesized map is obtained by coadding the ILC estimates at the various scales. The wavelet decomposition allows for estimating the dimension, m, of the Galactic signal subspace locally over the sky and over the angular scales, depending on the local conditions of contamination. In this work, the needlet bandpass windows (Fig. A.2) are defined in harmonic space from the difference of successive Gaussian beam transfer functions:
(A.29)where (A.30)and (A.31)with FWHM = [300′,120′,60′,45′,30′,15′,10′,7.5′,5′]. In this way we have (A.32)such that there is no effective filtering of any power at any scale in the final maps synthesized from the different needlet scales after component separation.
Fig. A.1 Local number of Galactic foreground degrees of freedom selected by the AIC criterion at one degree angular scale (top), 20′ scale (middle), and 5′ scale (bottom). The number of Galactic degrees of freedom decreases at high latitude and small angular scales. 
Fig. A.2 Needlet windows acting as bandpass filters in harmonic space (black lines), with the 5′ beam transfer function overplotted (red line). 
Appendix A.4: Model selection with the Akaike information criterion
In Remazeilles et al. (2011b), the effective number, m, of Galactic components in each needlet domain was estimated by rejecting the eigenvalues in Eq. (A.22) that are smaller than 1.25, i.e. for which the “noise” contributes to the observation by more than 80%. This criterion is somewhat arbitrary. In the present work, we propose instead to use a statistical criterion to discriminate between the “large” eigenvalues, tracing the Galactic signal, and the “noisy” eigenvalues (≈ 1) to be rejected; the effective rank of the covariance matrix of the Galactic signal is estimated by statistical model selection through the Akaike information criterion (Akaike 1974).
For a given dimension, or model, m, if we assume that the data, x, are independent and identically distributed according to the Gaussian distribution , with R(m) = R_{f}(m) + R_{N}, then the likelihood reads as (A.33)where n is the number of modes in the (needlet) domain considered. The loglikelihood can be written as (A.34)where is the KullbackLeibler divergence (Kullback 1968), measuring the spectral mismatch between the model covariance matrix, R(m), and the data covariance matrix, : (A.35)At this stage, it is interesting to note that the estimate of the Galactic covariance matrix, computed in Eq. (A.25), is nothing other than the maximum likelihood estimate, i.e. the minimizer of the KullbackLeibler divergence of Eq. (A.35), as in SMICA (Delabrouille et al. 2003; Cardoso et al. 2008). The proof is given in Sect. A.5.
In the region of the sky and the range of angular scales considered, we select the best rank value, m^{∗}, among the class of models, m, by minimizing the AIC (A.36)Through the penalty, 2 nm, the AIC makes a tradeoff between the goodness of fit and the complexity of the model. Let us denote the diagonalization of the transformed data covariance matrix, where (A.37)Since the maximum likelihood is reached when , we have (A.38)where μ_{i} are the (N_{ch}−m) eigenvalues of the transformed data covariance matrix, , collected in the matrix D_{N}, and (A.39)This convex function is minimum for μ = 1.
Therefore, the AIC criterion reduces to the simple analytical form: (A.40)and the dimension of the Galactic signal subspace is estimated by minimizing Eq. (A.40) in each region considered: (A.41)Minimizing the AIC criterion (Eq. (A.40)) in each needlet domain (i.e. in each region of the sky and each range of angular scales) allows for estimating the effective dimension of the Galactic signal subspace locally, given the level of contamination by noise, CIB, and CMB in this region. The multidimensional ILC then adapts the filtering to the effective local dimension of the Galactic signal. In this respect, the GNILC method goes beyond the SMICA method by relaxing any prior assumption on the number of Galactic components.
Figure A.1 shows the local dimension of the Galactic signal subspace over the sky, for three different ranges of angular scales, which we have estimated by minimizing the AIC criterion (Eq. (A.40)) on the Planck data. In practice, the dimension of the Galactic signal subspace (equivalently, the number of Galactic degrees of freedom) clearly depends on the local signaltonoise ratio, where the “noise” here includes the instrumental noise, the CIB, and the CMB signals. As expected, at high latitude and small angular scales (bottom panel in Fig. A.1) this number decreases because the CIB and the instrumental noise become dominant in the observations.
Appendix A.5: Minimization of the KullbackLeibler divergence
To close this appendix, we show that the GNILC estimates of the covariance matrix (Eq. (A.25)) and the mixing matrix (Eq. (A.27)) are the maximum likelihood solution that minimizes the KullbackLeibler divergence (Eq. (A.35)). For the sake of simplicity, we will assume that the data have been whitened so that the nuisance covariance matrix is represented by an identity matrix, I.
Let us consider a foreground model with a fixed number, m, of independent templates in the domain being considered. The foreground covariance matrix can therefore be modelled as a rank–m matrix (A.42)and the full data covariance matrix as (A.43)where Λ_{m} is an m × m diagonal matrix collecting the m nonnull eigenvalues of R_{f}(m) and A_{m} is an N_{ch} × m matrix collecting the m corresponding eigenvectors. In other words, the eigendecomposition of the foreground covariance matrix is given by
(A.44)The maximum likelihood estimate, Λ_{m}, minimizing the KullbackLeibler divergence, must satisfy for all (i,j) ∈ [1,m]^{2}(A.45)Therefore, the maximum likelihood estimate, { A_{m},Λ_{m} }, is solution of (A.46)where I is an m × m identity matrix. If we consider the eigenvalue decomposition of the full data covariance matrix as (A.47)where D_{S} is the m × m diagonal matrix collecting the m largest eigenvalues of , and U_{S} is an N_{ch} × m matrix collecting the m corresponding eigenvectors, then the maximum likelihood solution of Eq. (A.46) is given by (A.48)
Appendix A.6: Test on Planck simulations
We apply the GNILC algorithm to the public^{5}Planck full focal plane simulations (FFP8) from Planck Collaboration XII (2016) in order to validate the componentseparation method.
We still adopt as a prior for GNILC the Planck 2013 CIB bestfit power spectra, although the simulations of the CIB components in FFP8 do not follow the exact same statistics, since the FFP8 CIB maps are generated across frequencies from simulated dark matter shells in Planck Collaboration XII (2016). The imperfect agreement between the prior and the actual CIB power spectrum enables us to test the robustness of GNILC when the knowledge of the CIB power spectrum is not perfect.
In Fig. A.3 we compare the reconstructed GNILC dust map with the input FFP8 dust map within some small CIB fields of Planck Collaboration XXX (2014), where the dust is faint, namely NEP4 (part of the GHIGLS NEP field) and SP. The difference between the GNILC dust map and the input dust map is also shown. In Fig. A.4, we compare the reconstructed GNILC CIB map with the input FFP8 CIB map. The difference between the GNILC CIB and the input CIB is shown in the third column of Fig. A.4.
The T–T correlation between the GNILC output dust and the FFP8 input dust is plotted in Fig. A.5. The T–T correlation between the GNILC output CIB and the FFP8 input CIB is plotted in Fig. A.6. In the NEP4 field, the slope of the T–T scatter plots is 0.997 ± 0.006 for the dust and 0.981 ± 0.008 for the CIB. In the SP field, the slope of the T–T scatter plots is 0.997 ± 0.006 for the dust and 1.079 ± 0.008 for the CIB. In both fields, the Pearson correlation coefficient between the input and the GNILC reconstruction is found larger than 0.9 for the dust and larger than 0.8 for the CIB.
The successful reconstruction of dust and CIB shows that the GNILC method is robust to imperfect prior assumptions on the CIB angular power spectra. As shown in Olivari et al. (2016), it is not the exact morphology of the CIB power spectrum but the overall slope that matters for enabling GNILC to discriminate between dust and CIB. The prior power spectrum in GNILC is only intended to estimate the dimension of the dust and CIB component subspaces (Sect. A.3) not the amplitude of these components (Sect. A.1).
Fig. A.3 Thermal dust maps from the Planck FFP8 simulations. Input dust FFP8 (left), GNILC dust FFP8 (middle), difference (GNILC – input) (right). The fields are: NEP4 (top) and SP (bottom). 
Fig. A.4 CIB maps from the FFP8 simulations. Input CIB FFP8 (left), GNILC CIB FFP8 (middle), difference (GNILC – input) (right). The fields are: NEP4 (top) and SP (bottom). 
Fig. A.5 T–T scatter plots of the Planck FFP8 simulations between GNILC and the input thermal dust maps in the NEP4 field (left) and the SP field (right). 
Fig. A.6 T–T scatter plots of the Planck FFP8 simulations between GNILC and input CIB maps in the NEP4 field (left) and the SP field (right). 
All Tables
Mean and dispersion of the dust temperature and spectral index in different areas of the sky (full sky, high latitude, low latitude). Top: GNILC MBB fit. Middle: PR2 MBB fit a la dust model P13. Bottom: dust model P15 (Commander60′).
All Figures
Fig. 1 gnomonic projection of the sky centred at high latitude (l,b) = (90°,−80°). Top left: Planck 353GHz map. Top right: CMBremoved Planck 353GHz map. Middle left: dust model P13 at 353 GHz (MBB fit on CMBremoved Planck maps). Middle right: GNILC dust map at 353 GHz. Bottom left: dust model P15 at 545 GHz (Commander Bayesian fitting). Bottom right: GNILC dust map at 545 GHz. Maps at 353 GHz are shown at 5′ resolution, while maps at 545 GHz are smoothed to 7.5′ resolution. The GNILC dust maps have a nonuniform resolution (see Fig. 2) with 5′ resolution kept in regions of bright dust emission. In each image the local mean intensity has been subtracted for this comparison. 

In the text 
Fig. 2 Effective beam FWHM of the GNILC dust maps on the whole sky (top panel) and on a area of the sky centred at high latitude (l,b) = (90°,−80°) (bottom panel). The spatially varying beam FWHM is the same for all frequencies. GNILC preserves the 5′scale power of the thermal dust in the high signaltonoise regions of the sky. 

In the text 
Fig. 3 Fullsky map of the Galactic thermal dust emission: Planck 2013 (P13) thermal dust model at 353 GHz and 5′ resolution (top panel), suffering from CIB contamination at high latitudes, and the GNILC dust map (this work) at 353 GHz and 5′ resolution (bottom panel), for which the CIB is clearly filtered out at highlatitudes. A logarithmic colour scale is used here to highlight the lowintensity emission at high latitudes. The effective local beam of the GNILC dust maps is shown in Fig. 2. 

In the text 
Fig. 4 GNILC dust maps at 353 GHz (left panel), 545 GHz (middle panel), and 857 GHz (right panel) on a gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). GNILC filters out the CIB anisotropies while preserving the smallscale dust signal (see bottom panel of Fig. 2). In these images, the local mean intensity of each map has been subtracted. 

In the text 
Fig. 5 gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). Left: difference map (Planck 353 GHz −Planck CMB −GNILC dust) reveals the CIB anisotropies at 353 GHz. Middle: difference map (Planck 353 GHz – Planck CMB – dust model P13) revealing only the instrumental noise because the dust model P13, like the Planck observations at 353 GHz, still contains the CIB signal. Right: difference (dust model P13 – GNILC dust) revealing the amount of CIB contamination in the dust model P13 with respect to the GNILC dust map. In these images, the local mean intensity of each map has been subtracted. 

In the text 
Fig. 6 Angular power spectra of the various maps at 353 GHz (top panel), 545 GHz (middle panel), and 857 GHz (bottom panel), on a fraction of the sky, f_{sky} = 57%: Planck map (solid black line), dust model P13 (dotted blue line, Planck Collaboration XI 2014), GNILC dust map (long dashed red line), GNILC dust map corrected for the residual noise (solid green line), and GNILC residual map (Planck map −Planck CMB −GNILC dust, solid yellow line), which is compared to the Planck CIB bestfit power spectrum (dashdot purple line, Planck Collaboration XXX 2014) and the Planck noise power spectrum (dashed orange line). 

In the text 
Fig. 7 Fullsky thermal dust parameter maps: temperature (top row), spectral index (middle row), and map of the χ^{2} statistic of the fit (bottom row). Left panels: GNILC modified blackbody (MBB) fit. Right panels: PR2 modified blackbody (MBB) fit a la model P13. 

In the text 
Fig. 8 gnomonic projections of the dust temperature maps at high latitude, b = −80° (top panels) and low latitude, b = −20° (bottom panels). Left: GNILC MBB fit. Right: PR2 MBB fit according to model P13. 

In the text 
Fig. 9 gnomonic projections of the dust spectral index maps at high latitude, b = −80° (top panels) and low latitude, b = −20° (bottom panels). Left: GNILC MBB fit. Right: PR2 MBB fit according to model P13. 

In the text 
Fig. 10 Highlatitude area of the sky with f_{sky} = 21% (top) and lowlatitude area of the sky with f_{sky} = 20% (bottom) that are considered in Fig. 11 and Table 1. 

In the text 
Fig. 11 Normalized histograms of T_{dust} and β_{dust} at 5′ resolution for the GNILC MBB fit (red contours) and the PR2 MBB fit a la model P13 (green contours). The normalized histograms for the dust model P15 (Commander fit at 60′ resolution) are overplotted (blue contours). The histograms are computed from the subset of pixels corresponding to either the highlatitude area in the sky with f_{sky} = 21% (upper panels), the lowlatitude area in the sky with f_{sky} = 20% (middle panels), or the whole sky (lower panels). Due to CIB contamination at highlatitude, the PR2 MBB fits show larger dispersion than the GNILC MBB fits in the distributions of T_{dust} and β_{dust}. 

In the text 
Fig. 12 GNILC CIB maps for a large fraction of the sky at 353 GHz (top), 545 GHz (middle), and 857 GHz (bottom). Apart from thermal dust reconstruction, the GNILC componentseparation method gives access to CIB anisotropies over 57% of the sky. The left panels show the GNILC CIB maps at full resolution while the right panels show the CIB maps smoothed to one degree resolution. 

In the text 
Fig. 13 GNILC CIB maps at 353 GHz (left), 545 GHz (middle), and 857 GHz (right) on a gnomonic projection of the sky centred at high latitude, (l,b) = (90°,−80°). The partial spatial correlation of the CIB anisotropies between pairs of frequencies is clearly visible. In this figure, the local mean intensity of each map has been subtracted. 

In the text 
Fig. 14 Maps of cosmic infrared background (CIB) anisotropies in GHIGLS fields (Martin et al. 2015). Left column: CIB 2013 (Planck Collaboration XXX 2014). Middle column: CIB GNILC. Right column: difference (CIB 2013 − GNILC). From top to bottom row, the GHIGLS fields are: Boötes; MC; N1; and SP. The size and the location in the sky of each field are defined in Table 1 of Martin et al. (2015). 

In the text 
Fig. 15 T–T scatter plot between CIB GNILC and CIB 2013. The fields are Boötes, MC, N1, and SP. 

In the text 
Fig. 16 T–T scatter plot between the H i map and the GNILC CIB map at 353 GHz over 57% of the sky (see Fig. 12). 

In the text 
Fig. 17 T–T scatter plot at high Galactic latitude between the H i map and the dust optical depth maps at 353 GHz: dust model P13 (green), PR2 MBB fit (black), GNILC MBB fit (blue). While the slope of the correlation with the H i map is consistent for all the dust optical depth maps, the scatter is smallest for the GNILC dust optical depth map because of the removal of the CIB temperature anisotropies. 

In the text 
Fig. 18 Polar orthographic projections of the E(B−V) maps at north (left) and south (right) Galactic poles in the PanSTARRS 1 survey area (Onaka et al. 2008). Top: E(B−V) map from Green et al. (2015). Bottom: GNILC E(B−V) map. 

In the text 
Fig. 19 T–T scatter plot between the ratio E(B−V) /N_{H} and the gas column density N_{H} for the twodimensional projection of the E(B−V) map of Green et al. (2015) (red diamonds) and the GNILC E(B−V) map (blue diamonds) in the highlatitude region of the sky defined in Fig. 10. Each point is the average of E(B−V) /N_{H} values in a bin of N_{H}. The bin size varies such that there is always the same number of samples per bin. 

In the text 
Fig. A.1 Local number of Galactic foreground degrees of freedom selected by the AIC criterion at one degree angular scale (top), 20′ scale (middle), and 5′ scale (bottom). The number of Galactic degrees of freedom decreases at high latitude and small angular scales. 

In the text 
Fig. A.2 Needlet windows acting as bandpass filters in harmonic space (black lines), with the 5′ beam transfer function overplotted (red line). 

In the text 
Fig. A.3 Thermal dust maps from the Planck FFP8 simulations. Input dust FFP8 (left), GNILC dust FFP8 (middle), difference (GNILC – input) (right). The fields are: NEP4 (top) and SP (bottom). 

In the text 
Fig. A.4 CIB maps from the FFP8 simulations. Input CIB FFP8 (left), GNILC CIB FFP8 (middle), difference (GNILC – input) (right). The fields are: NEP4 (top) and SP (bottom). 

In the text 
Fig. A.5 T–T scatter plots of the Planck FFP8 simulations between GNILC and the input thermal dust maps in the NEP4 field (left) and the SP field (right). 

In the text 
Fig. A.6 T–T scatter plots of the Planck FFP8 simulations between GNILC and input CIB maps in the NEP4 field (left) and the SP field (right). 

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.