© ESO, 2016

1. Introduction

This paper, one of a set associated with the 2015 release of data from the Planck¹ mission (Planck Collaboration I 2016), presents a coherent astrophysical model of the microwave sky in both temperature and polarization, as derived from the most recent Planck observations. For temperature, the analysis also incorporates the 9-yr WMAP observations (Bennett et al. 2013) and a 408 MHz survey (Haslam et al. 1982), allowing the separation of synchrotron, free-free, and spinning dust emission.

In March 2013, the Planck Consortium released its first temperature measurements of the microwave sky, summarized in terms of nine frequency maps between 30 and 857 GHz (Planck Collaboration I 2014). The richness of these data has enabled great progress in our understanding of the astrophysical composition of the microwave sky. The current Planck data release presents additionally high-sensitivity, full-sky maps of the polarized microwave sky, offering a fresh view on both cosmological and astrophysical phenomena.

With increased data volume and quality comes both greater scientific potential and more stringent requirements on model complexity and sophistication. The current Planck data release is more ambitious than the 2013 release in terms of component separation efforts, accounting for more astrophysical effects and components. In this round, three related papers summarize the Planck 2015 component separation products and approaches. First, cosmic microwave background (CMB) reconstruction and extraction are discussed in Planck Collaboration IX (2016). Second, this paper presents the diffuse astrophysical foreground products derived from the 2015 Planck observations, both in temperature and polarization. Third, Planck Collaboration XXV (2016) discusses the scientific interpretation of the new low-frequency Planck foreground products.

The main goal of the current paper is to establish a single, internally coherent and global parametric model of the microwave sky, simultaneously accounting for all significant diffuse astrophysical components and relevant instrumental effects using the Bayesian Commander analysis framework (Eriksen et al. 2004, 2006, 2008). As such, our discussion does not focus on any single emission component, but rather emphasize the global picture. In the 2013 data release, the same framework was applied to the Planck temperature measurements for frequencies between 30 and 353 GHz, considering only angular scales larger than 40′ full-width half-maximum (FWHM). This resulted in low-resolution CMB, CO, and thermal dust emission maps, as well as a single low-frequency foreground component combining contributions from synchrotron, free-free, and spinning dust emission (Planck Collaboration XII 2014). Here we extend that analysis in multiple directions. First, instead of 15.5 months of temperature data, the new analysis includes the full Planck mission data, 50 months of Low Frequency Instrument (LFI) and 29 months of High Frequency Instrument (HFI) data, in both temperature and polarization. Second, we now also include the 9-yr WMAP observations between 23 and 94 GHz and a 408 MHz survey map, providing enough frequency constraints to decompose the low-frequency foregrounds into separate synchrotron, free-free, and spinning dust components. Third, we now include the Planck 545 and 857 GHz frequency bands, allowing us to constrain the thermal dust temperature and emissivity index with greater precision, thereby reducing degeneracies between CMB, CO, and free-free emission. At the same time, we find that the calibration and bandpass measurements of these two channels represent two of the most important sources of systematic uncertainty in the analysis. Fourth, the present analysis implements a multi-resolution strategy to provide component maps at high angular resolution. Specifically, the CMB is recovered with angular resolution 5′ FWHM (Planck Collaboration IX 2016), thermal dust emission and CO J = 2 → 1 lines are recovered at 7.́5 FWHM, and synchrotron, free-free, and spinning dust are recovered at 1deg FWHM. The resulting parameter fits define the Planck 2015 baseline astrophysical model in temperature and polarization. We emphasize, however, that these models are not unique, but instead represent minimal physically well-motivated models that are able to reproduce the current data.

As in the 2013 data release, the CMB solutions derived, using this Bayesian approach, form the basis of the Planck 2015 CMB temperature likelihood on large angular scales. This is described in detail in Planck Collaboration XI (2016), which also presents a detailed characterization of the low-multipole CMB angular power spectrum. The low-frequency astrophysical model presented here is used as input for the temperature-to-polarization bandpass mismatch corrections for the LFI polarization maps (Planck Collaboration II 2016).

The paper is organized as follows. Section 2 gives an overview of the computational framework implemented in the Commander code. Section 3 describes the data selection and processing. Section 4, gives an overview of the relevant astrophysical components and systematic effects. Sections 5 and 6 give the main temperature and polarization products. We summarize in Sect. 7.

2. Algorithms

2.1. Data, posterior distribution and priors

Most of the results derived in this paper are established within a standard Bayesian analysis framework, as implemented in the Commander code, in which an explicit parametric model, s(θ), is fitted to a set of observations, d, either by maximizing or mapping out the corresponding posterior distribution, (1)Here θ denotes some general set of free parameters in the model, ℒ(θ) = P(d | θ) is the likelihood, and P(θ) denotes a set of priors on θ. The evidence, P(d), is a constant with respect to the parameter set, and is neglected in the following.

The data are defined by a set of pixelized frequency-channel sky maps, d = { d_ν }, comprising the three Stokes parameters I, Q and U. In this paper, however, we analyse temperature and polarization separately; therefore the data vector comprises either I or { Q,U }.

We start by assuming that the data at a given frequency ν may be described as a linear sum of signal s_ν and noise n_ν, (2)where n_ν is assumed to be Gaussian-distributed with a known covariance matrix N_ν. For the signal, we adopt the following parametric expression: where a_i is an amplitude map for component i at a given reference frequency, β_i is a general set of spectral parameters for the same component, g_ν is a multiplicative calibration factor for frequency ν, Δ_ν is a linear shift in the bandpass central frequency, and m_ν is a set of template correction amplitudes, such as monopole, dipole, or zodiacal light corrections for temperature, or calibration leakage templates for polarization. The corresponding spatially fixed templates are organized column-wise in a template matrix T_ν. The mixing matrix, , accounts for the effect of spectral changes as a function of frequency for component i, parametrized by β_i, as well as bandpass integration effects and unit conversions. For numerical stability, all internal calculations are performed in units of brightness temperature, and a is therefore naturally defined in the same units at some specified reference frequency.

The posterior distribution takes the usual form, (5)where we have included the CMB power spectrum, C_ℓ, and also implicitly adopted uniform priors on g_ν, Δ_ν, m_ν, and C_ℓ. Because the noise is assumed to be Gaussian and independent between frequency channels, the likelihood reads (6)Likewise, we further assume the CMB signal to be Gaussian distributed with a covariance matrix, S(C_ℓ), given by the power spectrum, and the corresponding CMB prior factor therefore reads (7)

The only undefined factors in the posterior are the amplitude and spectral parameter priors, P(aⁱ) and P(βⁱ). These represent the most difficult problem to handle from a conceptual point of view, since the prior is to some extent a matter of personal preference. However, we adopt the following general practices in this paper. First, for low-resolution analyses that include fitting of template amplitudes (e.g., monopoles and dipoles), we always impose a strict positivity prior, i.e., a_i> 0, on all signal amplitudes except the CMB. Without such a prior, there are large degeneracies between the zero-points of the amplitude maps and the individual template amplitudes. Second, for the high angular resolution analysis, we fix the template amplitudes at the low-resolution values and disable the positivity prior, in order to avoid noise bias. Third, to further break degeneracies, we adopt fiducial values for the monopole, dipole, and calibration factors for a few selected channels, effectively imposing a set of external priors from CMB dipole measurements and H i cross-correlation (Planck Collaboration VIII 2014) to anchor the full solution. Fourth, for the spectral parameters, we adopt Gaussian priors with means and variances informed by the high signal-to-noise values observed in the Galactic plane, which for all practical purposes are independent of the adopted priors. Intuitively, we demand that a map of the spectral parameter in question should not be much different in the data-dominated and the prior-dominated regions of the sky. Fourth, one of the components in the temperature model is free-free emission, which has two free parameters, namely the effective emission measure, EM, and the electron temperature, T_e. The latter of these is very poorly constrained with the current data set except in the central Galactic plane, and we therefore adopt a smoothness prior on this paper to increase the effective signal-to-noise ratio, demanding that is must be smooth on 2deg FWHM scales. This in turn has a large computational cost by making the overall foreground parameter estimation process non-local, and T_e is therefore only varied in fast maximum-likelihood searches, not in expensive sampling analyses. Its effect on other parameters is, however, minimal, precisely because of its low signal-to-noise ratio. Finally, in addition to these informative priors, we adopt a Jeffreys prior for the spectral parameters in order to suppress prior volume effects (Jeffreys 1946; Eriksen et al. 2008; Dunkley et al. 2009).

2.2. Gibbs sampling and posterior maximization

As described above, the posterior distribution contains many millions of free (both non-Gaussian and strongly correlated) parameters for Planck – ranging from 11 million in the following low-resolution analysis to 200 million in the corresponding high-resolution analysis – and mapping out this distribution poses a significant computational problem. Indeed, no direct sampling algorithm exists for the full distribution, and the only computationally efficient solution currently known is that of Gibbs sampling, a well-known textbook algorithm in modern statistical analysis (e.g., Gelman et al. 2003). The underlying idea of this method is that samples from a complicated multivariate distribution may be drawn by iteratively sampling over the corresponding conditional distributions, which usually have much simpler, and often analytic, sampling algorithms. This framework was originally introduced to the CMB analysis field by Jewell et al. (2004) and Wandelt et al. (2004), and subsequently developed into a fully functional computer code called Commander by Eriksen et al. (2004, 2008).

For the problem in question in this paper, this algorithm may be schematically translated into an explicit set of sampling steps through the following Gibbs chain: Here “←” denotes drawing a sample from the distribution on the right-hand side. After some burn-in period, the theory of Gibbs sampling guarantees that the joint set of parameters is indeed drawn from the correct joint distribution. For a full description of the various steps in the algorithm, see Eriksen et al. (2008).

While no fully functional alternatives to Gibbs sampling have been established for this full joint distribution to date, Gibbs sampling alone by no means solves all computational problems. In particular, this algorithm is notorious for its slow convergence for nearly degenerate parameters, since it by construction only moves through parameter space parallel to coordinate axes. For this reason, we implement an additional posterior maximization phase, in which we search directly for the posterior maximum point rather than attempt to sample from the full distribution. The resulting solution may then serve either as a final product in its own right, by virtue of being a maximum-posterior estimate, or as the starting position for a regular Gibbs sampling analysis. The crucial point, though, is that special-purpose nonlinear search algorithms can be used in this phase, moving in arbitrary directions through parameter space, and individual optimization combinations may be introduced to jointly probe directions with particularly strong degeneracies. Perhaps the single most important example in this respect is the parameter combination between component amplitudes, detector calibrations, and bandpass uncertainties, { a_i,g_ν,Δ_ν }, all three of which essentially correspond to scaling parameters. However, since both g_ν and a_i are conditionally linear parameters, and only Δ_ν is truly nonlinear, it is possible to solve analytically for g_ν or a_i, conditioning on any given fixed value of Δ_ν. Consequently, one can set up a non-linear Powell-type search (Press et al. 2002) for Δ_ν, in which the optimal values of either g_ν or a_i are quickly computed at each iteration in the search. A second example is the electron temperature discussed above, for which non-local optimization is feasible, whereas a full-blown sampling algorithm is too expensive to converge robustly. In this situation, fixing the parameter at its maximum-posterior value is vastly preferable compared to adding an unconverged degree of freedom in the full sampler.

Even with this optimization phase, however, there is always an inherent danger of the algorithm being trapped in a local posterior maximum. Indeed, with a distribution involving millions of highly correlated parameters, it is exceedingly difficult to prove that the derived solution is the true global posterior maximum. As a partial solution to this problem, we initialize the search using different starting positions, and carefully monitor the convergence properties of the chains.

3. Data selection and processing

The primary data used in this paper are the 2015 Planck temperature and polarization sky maps (Planck Collaboration VI 2016; Planck Collaboration VIII 2016). For the temperature analysis we additionally include the 9-yr WMAP observations² (Bennett et al. 2013) and a full-sky 408 MHz survey map (Haslam et al. 1982), with the goal of individually resolving synchrotron, free-free, and spinning dust emission. For WMAP, we adopt the beam-symmetrized frequency maps for the foreground-dominated K- and Ka-bands, to mitigate beam artifacts around compact sources, but we use the standard maps for the CMB-dominated Q-, V-, and W-bands, because of their more accurate noise description. At the lowest frequency, we adopt the destriped version of the 408 MHz survey map recently published by Remazeilles et al. (2015).

In order to maximize our leverage with respect to bandpass measurement uncertainties and line emission mechanisms, we employ individual detector and detector set (“ds”) maps for all Planck frequencies between 70 and 857 GHz, and differencing assembly (“DA”) maps for WMAP. However, the polarization analysis employs frequency maps in order to maximize signal-to-noise ratio and to minimize correlated noise from destriper mapmaking uncertainties. Intensity-to-polarization leakage from bandpass mismatch between detectors is suppressed through the use of precomputed leakage templates (Planck Collaboration II 2016; Planck Collaboration VIII 2016). For LFI, these templates are based on a preliminary version of the foreground products presented in this paper. The full set of clean channels used in this analysis is summarized in Table 1 in terms of centre frequencies, resolution, and noise levels. In total, 32 individual detector and detector set maps³ and frequency maps are included in the temperature analysis, and seven frequency maps in the polarization analysis.

For the Planck HFI channels, a model of zodiacal light emission is subtracted from the time-ordered data prior to mapmaking (Planck Collaboration VIII 2016). In addition, in this paper we apply a small correction to the low-frequency LFI and WMAP channels by scaling the effective HFI 100 GHz zodiacal light correction map (i.e., uncorrected minus corrected map) to each frequency according to a power law fitted to frequencies between 100 and 353 GHz (Planck Collaboration XIV 2014), as illustrated in Fig. 1; the actual template amplitudes relative to the 100 GHz correction map (in thermodynamic units) are listed in Table 2. Although the magnitude of this correction is small, with a maximum amplitude of 2μK in the 70 GHz map, applying no correction at all below 100 GHz results in a visually noticeable bias in the derived CO J = 1 → 0 map at high Galactic latitudes, in the characteristic form of the zodiacal light. Extending the zodiacal light model to low frequencies efficiently eliminates this structure.

In our 2013 release, colour corrections and unit conversions for all Planck channels were based on individual bandpass profiles as measured on the ground before launch (Planck Collaboration V 2014; Planck Collaboration IX 2014). However, as discussed in detail in Sects. 2, 4.3, and 5, during the component separation process we find that systematic uncertainties in the nominal bandpasses induce significant residuals between data and model, and it is necessary to fit for these bandpass uncertainties in order to obtain statistically acceptable fits. For WMAP, we adopt the nominal bandpasses⁴ for the first DA within each frequency band, and fit for the remaining DA bandpasses within each frequency. For the 408 MHz survey, we adopt a delta function response at the nominal frequency. The unit conversion factors between thermodynamic and brightness temperatures and flux density per area are tabulated for both nominal and fitted bandpass profiles in Table 3.

Some Planck detector maps are affected more strongly by systematic errors than others (Planck Collaboration VIII 2016), and in the following temperature analysis we exclude the worst channels in order not to compromise the overall solution. Out of a total of 31 potential Planck detector and detector set maps, 10 are removed from further analysis, while the remaining 21 are listed in Tables 1 and 3. The 10 removed maps are shown in Fig. 2 in the form of a (d_ν−s_ν) residual map, where s_ν is a signal model based on the global and spectral parameters derived in Sect. 5, but with amplitudes re-fitted to all 42 channels. The top panel shows the corresponding χ² map, defined as (14)Including any one of these ten maps in the full joint analysis increases the χ² of the total fit far beyond what is allowed by random statistical fluctuations.

Several different classes of systematic effects may be seen in these plots. Starting with the 217-ds1 map, we see a large-scale red-blue pattern, aligned with the Planck scanning strategy and crossing through the Ecliptic poles (Planck Collaboration I 2014). The same spatial pattern is seen in the 217-ds2, 353-ds1, 353-2, 353-7, 353-8, 857-3, and 857-4 maps as well, with varying amplitudes and signs. These are due to low-amplitude destriping errors induced in the mapmaking process (Planck Collaboration VIII 2016), and only become visually apparent after component separation removes the dominant Galactic signal. These residuals may therefore, at least in principle, be suppressed by iterating between mapmaking and component separation, essentially performing a joint mapmaking/component separation χ² fit for the destriping offsets.

A second example of residual systematics is seen in the 857-1 and 857-4 maps, and to a lesser extent in the 545-1 map, in the form of sharp features at high Galactic latitudes. These correspond to far sidelobe (FSL) contamination at the 0.05 MJy sr^-1 level. Third, and somewhat more subtly, one can see the effect of the very-long time constants (VLTC) discussed in Planck Collaboration VIII (2016), particularly south of the Galactic plane in 353-2, 353-7, and 545-1. In each of these maps, the Galactic plane appears to have been smeared along the scanning direction. As in the case of destriping errors, combining mapmaking and component separation should prove very powerful in suppressing both FSL and VLTC errors in future analyses.

While several channels are omitted from the analysis, it is important to note that the Planck HFI temperature observations are strongly signal-dominated at all angular scales above 7′ FWHM. Our data cuts therefore have very little effect on the full error budget, which is dominated by modelling and systematic errors for most parameters. The main cost of these cuts comes in the form of reduced internal redundancy. In particular, having access to only one clean 857 GHz channel limits our ability to determine its bandpass and calibration coefficients. This situation will improve in the next Planck data release, when the remaining three 857 GHz channels are better cleaned at the time domain level.

The current analysis is carried out in a number of stages according to angular resolution, in which higher-resolution stages implement a simpler foreground model than lower-resolution stages. For temperature, each full-resolution map is downgraded from its native resolution to 1deg (all channels), 40′ (all Planck channels), and (HFI channels above 100 GHz) by deconvolving the intrinsic instrumental beam profile and convolving with a Gaussian beam of the appropriate size before repixelizing at HEALPix⁵ resolutions N_side = 256, 256, and 2048, respectively (Górski et al. 2005). For polarization, the corresponding smoothing scales are 40′ and 10′ FWHM, pixelized at N_side = 256 and 1024, respectively.

The instrumental noise is assumed to be spatially uncorrelated and Gaussian for all channels, with a spatially-varying rms given by the scanning strategy of the experiment and the instantaneous sensitivity of each detector. The low-resolution noise rms map is found by convolving the high-resolution rms with the appropriate smoothing kernel, properly accounting for its matrix-like nature, and retaining only the diagonal element of the resulting covariance matrix. For the 545 and 857 GHz channels, we additionally add 0.01 MJy sr^-1 of uniform white noise, to prevent known residual far sidelobe and destriping contamination from propagating to lower frequencies through the thermal dust temperature and spectral index and contaminating both the CMB and CO solutions. Similarly, for the 408 MHz channel we add regularization noise, equal to 1% of the amplitude of the map, to account for low-level residuals not captured by the white noise model described above.

Finally, two different processing masks are employed in the temperature analysis (PM99.6 and PM61), removing 0.4% and 39% of the sky, respectively, as shown in Fig. 3. PM99.6 is generated by thresholding the (smoothed) χ² map of a preliminary analysis⁶ at 10⁴, removing only the very brightest outliers in the data set; this mask is used for bandpass estimation and gain calibration of the 545-4 and 857-2 channels. PM61 is generated as the product of the χ² map resulting from an analysis without bandpass corrections thresholded at the 5σ level, and the 9-yr WMAP point source mask. This mask is used for calibration estimation of the CMB frequencies, and for monopole and dipole estimation.

4. Model survey

Before presenting the results of our analysis, it is useful to review the various model components that are relevant for this work, with the goal of building intuition concerning important degeneracies and residuals that may be observed in various goodness-of-fit statistics. We first review each of the astrophysical sky signal components for both temperature and polarization, as summarized in Table 4 and Fig. 4, and then each of the most important instrumental effects.

4.1. Sky components in intensity

CMB – The CMB is given by a perfect blackbody with only a single spectral parameter, namely the CMB temperature. We adopt a mean value of T_cmb = 2.7255 ± 0.0006 K (Fixsen 2009), and notice that the uncertainty in this value is sufficiently small to justify its use as a delta function prior. In the current paper, we neglect both Rayleigh scattering and higher-order relativistic effects (Planck Collaboration XXVII 2014; Lewis 2013); these could be accounted for in future work. The resulting CMB brightness temperature is plotted as a dashed line in each panel of Fig. 4, with an amplitude of 70μK, corresponding to the CMB rms at 1° FWHM resolution, providing a useful consistent visual reference for other components.

Synchrotron – Diffuse synchrotron emission is generated by relativistic cosmic-ray electrons spiraling in the Galactic magnetic field. This radiation may be highly polarized, with a maximum polarization fraction of 0.75. Both theoretical models and observations suggest that the synchrotron spectrum is well approximated by a power law with an index of β_s ≈ −3 at frequencies above 20 GHz, but with significant flattening at low frequencies. In this paper we adopt a fixed spectral template for the synchrotron spectrum. We use a spectrum extracted from the GALPROP z10LMPD_SUNfE synchrotron model from Orlando & Strong (2013), as described in Planck Collaboration XXV (2016). For brevity we refer to this as the GALPROP model from now on. A description of the GALPROP⁷⁸ code can be found in Moskalenko & Strong (1998), Strong et al. (2007), Orlando & Strong (2013) and references therein. We allow this spectrum to be rigidly shifted in log ν–log S space with a single global frequency shift parameter, α, for the full sky; see Table 4 for the explicit definition. The main effect of such translation, however, is to modify the synchrotron amplitude at 408 MHz, leaving β_s at frequencies above 20 GHz essentially constant and equal to −3.11 for any realistic shift parameter; see Fig. 4. Thus, with the current data the synchrotron amplitude is determined almost exclusively by the 408 MHz survey, while the frequency spectrum is determined by the GALPROP model, with the only free spectral parameter being the relative normalization between the 408 MHz and higher frequency channels. Several alternative models (straight, broken, or logarithmically-curved power laws) were also considered in the preparation of this paper, but we found no statistically robust evidence in favour of any of them over the GALPROP spectrum. With only one very low frequency channel, the current data set contains very little information about the synchrotron spectral index. To constrain the synchrotron component, and break remaining degeneracies between the synchrotron, free-free, and spinning dust components, data from up-coming low-frequency experiments such as S-PASS (2.3 GHz; Carretti et al. 2009), C-BASS (5 GHz; King et al. 2010, 2014), and QUIJOTE (10–40 GHz; Rubiño-Martín et al. 2012) are critically important. For now, we prefer the GALPROP model over various power-law variations, simply because it provides an acceptable fit to the data with essentially no free additional parameters, and that it is based on a well-defined physical model. For further discussion of synchrotron emission in the Planck observations, see Planck Collaboration XXV (2016).

Free-free – Free-free emission, or bremsstrahlung, is radiation from electron-ion collisions, and consequently has a frequency spectrum that is well-constrained by physical considerations (Dickinson et al. 2003). We adopt the recent two-parameter model of Draine (2011), in which the free parameters are the (effective beam-averaged) emission measure, EM (i.e., the integrated squared electron density along a line of sight), and the electron temperature, T_e. As seen in Fig. 4, the spectrum is close to a power law at frequencies above 1 GHz, but exhibits a sharp break at lower frequencies, where the medium becomes optically thick and the brightness temperature becomes equal to the electron temperature. Intuitively, over our frequency range, the EM determines the amplitude of the free-free emission, scaling the spectrum up or down, while the electron temperature changes the effective power law index very slightly, from β = −2.13 for T_e = 500 K to β = −2.15 for T_e = 20 000 K. At low frequencies, free-free emission is significantly degenerate with synchrotron emission, but because its power-law index is flatter than for synchrotron, it extends into the high signal-to-noise HFI channels. As a result, it is in fact possible to measure the electron temperature for particularly high signal-to-noise regions of the sky when imposing the smoothness prior discussed in Sect. 2.1. The main difficulty regarding free-free emission lies in a four-way degeneracy between CMB, free-free, spinning dust, and CO emission.

Thermal dust – At frequencies above 100 GHz, the dominant radiation mechanism is thermal dust emission Planck Collaboration XI (2014; Planck Collaboration XIV 2014; Planck Collaboration XVII 2014; Planck Collaboration Int. XIX 2015). The characteristic frequency is determined by the temperature of the dust grains, and therefore varies with dust population and environment. Empirically, thermal dust emission may be accurately described across the Planck frequencies as a modified blackbody with a free emissivity index, β_d, and dust temperature, T_d, per pixel, often referred to as a greybody, (Planck Collaboration XII 2014; Planck Collaboration Int. XVII 2014). At frequencies above 857 GHz, extending into the COBE-DIRBE frequencies (Hauser et al. 1998), the dust physics quickly becomes far more complicated, and the instrumental systematics more difficult to contain, and we therefore restrict our current analysis to frequencies up to 857 GHz. Correspondingly, we emphasize that the model derived here is only expected to be accurate at frequencies up to 857 GHz, and its main application is extrapolation to low, not high, frequencies. However, alternative parametrizations of thermal dust have been suggested in the literature, for instance in terms of a two-component grey-body model (e.g., Finkbeiner et al. 1999; Meisner & Finkbeiner 2015), in which each of the two components has a spatially constant emissivity and temperature. Similar to the alternative synchrotron models mentioned above, we also considered such models in the preparation of this paper. Without higher-frequency observations we found that the data are not able to discriminate between the various scenarios at any statistically significant level. For simplicity, we therefore adopt the one-component greybody model with a varying spectral index as a phenomenological parametric thermal dust model, noting that it provides a highly efficient representation of the true sky for frequencies up to 857 GHz, with no additional free global parameters. Finally, we point out that cosmic infrared background fluctuations (CIB) is interpreted as thermal dust emission in our model, by virtue of having a frequency spectrum similar to a greybody component (Planck Collaboration XXX 2014). However, because CIB is relatively weak compared to the Galactic signal, this is only a significant issue in low-foreground regions of the sky.

Spinning dust – Dust grains not only vibrate, but they can also rotate. If they have a non-zero electric dipole moment, this rotation will necessarily lead to microwave emission, as first demonstrated by Erickson (1957) and, in detail, by Draine & Lazarian (1998). The first hints of such radiation in real data were found by cross-correlating CMB data with far-IR observations (Kogut et al. 1996; de Oliveira-Costa et al. 1997; Netterfield et al. 1997), but at the time this correlation was thought to be dominated by free-free emission, coupled to dust due to star formation. The first identification of a new anomalous component was made by Leitch et al. (1997), and the observed excess emission was simply referred to as “anomalous microwave emission”, in order not to make premature conclusions regarding its physical nature. However, recent observational progress has made the physical association much more firm (e.g., Davies et al. 2006; Planck Collaboration XX 2011; Gold et al. 2011; Planck Collaboration XXI 2011; Bennett et al. 2013; Planck Collaboration Int. XV 2014, and references theirin) and we refer to this component simply as spinning dust. We adopt a spinning dust spectral template as our spectral model, as evaluated for a cold neutral medium with SpDust2 (Ali-Haïmoud et al. 2009; Ali-Haïmoud 2010; Silsbee et al. 2011). To fit this template to the data, we introduce a free peak frequency, ν_p, and allow rigid translations the log ν–log S space. As explicitly demonstrated by Bennett et al. (2013), this simple two-parameter model can accommodate a large number of spinning dust model variations to high precision.

To facilitate comparison with previously reported results, we report the peak frequency for a spinning dust spectrum defined in flux density units, and note that the conversion between brightness temperature and flux density is ∝λ^-2; see Table 4. The exact value of ν_p is highly uncertain, and depends on the physical properties of the dust grains, but is typically found to lie between 20 and 30 GHz from dedicated observations of individual objects (Planck Collaboration XX 2011), in some objects reaching as high as ≈50 GHz (Planck Collaboration Int. XV 2014). The exact shape of the spectrum is also uncertain, but we have been unable to fit this component with a single spdust component, always finding significant residuals with dust-like morphology at 20–30 GHz when forcing a high prior on ν_p, or at 40–70 GHz when forcing a low prior. For now, therefore, we adopt a two-component spinning dust model in which one component has a free peak frequency, ν_p, per pixel and the other component has a free, but spatially constant, ν_p2, for a total of three free spinning dust parameters per pixel and one global parameter. We emphasize that this is a purely phenomenological model, and we attach no physical reality to it. The second component is in effect only a computationally convenient parameterization of the width of the spinning dust peak (see Fig. 4). We therefore show only the sum of the two spinning dust components in the following, evaluated at 30 GHz. However, the set of released foreground data products includes the individual component parameters, allowing external uses to reproduce all results.

As listed in Table 4, we adopt a Gaussian prior of N(19 GHz ± 3 GHz) on the peak frequency of the primary spinning dust component, while no informative prior is enforced on the spatially uniform secondary component. At first sight, the former prior might look discrepant with respect to other targeted studies of spinning dust, which often report best-fit peak frequencies between 25 and 30 GHz. However, when adding the two spinning dust components, one with ν_p ≈ 20 GHz and another with ν_p ≈ 33 GHz, the total is indeed a component with a peak frequency around 25 GHz.

Finally, we emphasize that the spinning dust component exhibits significant correlations with both synchrotron and free-free emission, and its individual parameters are sensitive to both instrumental bandpasses and calibration. Thus, the reported parametric best-fit values are associated with large correlated systematic uncertainties, and additional low-frequency observations are required to make the spinning dust model robust. As already mentioned, low-frequency observations from experiments like S-PASS, C-BASS, and QUIJOTE are critically needed to break these degeneracies.

CO and 94/100 GHz line emission – As reported in Planck Collaboration XIII (2014), line emission from carbon monoxide (CO) is strongly detected in the Planck 100, 217, and 353 GHz frequency channels, and two sets of individual CO J = 1 → 0, 2 → 1, and 3 → 2 line maps (so-called “Type-1” and “Type-2” maps, produced by the MILCA and Commander algorithms, respectively) as well as a frequency co-added line map (“Type-3”, also produced with Commander) were provided in the 2013 Planck data release. In Sect. 5.5.3 we present updated CO maps, based on the full Planck mission data set, plus a new general line emission map, which captures a combination of emission lines that are detected with the HFI 100 GHz and WMAP W bands. An important contributor to this map is the HCN line at 88.6 GHz, providing about 23% (63%) of the 100 GHz (W-band) amplitude towards the Galactic circumnuclear disk (GCD) and Sgr A^∗ (Takekawa et al. 2014). Several other lines (CN, HCO+, CS, etc.) contribute at a level of 5–10% each. In addition, since we account for neither velocity effects nor detailed cloud physics (opacity, local thermal equilibrium state, etc.), there is also a non-negligible amount of CO leakage in this new map. We therefore refer to the new component simply as “94/100 GHz line emission”.

Generally speaking, the fact allowing us to separate line emission from other diffuse components is that the bandpass filters of individual detectors have different transmission levels at the line frequency. Therefore, while components with continuous spectra very nearly vanish in detector difference maps within single frequencies, line emission does not. For this reason, we employ individual detector maps in the 2015 temperature analysis, as opposed to co-added frequency maps in the corresponding 2013 analysis. However, for the highest resolution maps, for which neither LFI nor WMAP are able to provide useful information, there is a strong degeneracy between CMB, free-free, and CO J = 1 → 0 emission. The high-resolution CO maps are therefore produced with dedicated and special-purpose methods, as described in Sect. 5.5.3. Furthermore, the CO J = 3 → 2 emission map is significantly degenerate with the brighter thermal dust component, and consequently subject to large systematic uncertainties.

We describe the line emission maps parametrically in terms of an amplitude map at the line frequency, a(p), normalized relative to one specific detector map, and with a rigid frequency scaling given by the product of a spatially constant line ratio, h_ij, where i corresponds to a spectral line index, and j denotes detector index. F(ν) denotes some unit conversion factor, converting for example between μK_RJ and K_RJ km s^-1 (see Table 4 for the exact mathematical definition). The same formalism applies whether h_ij refers to line ratios between frequencies or only to bandpass ratios within frequencies. For example, in 2013 this formalism was used to construct the so-called CO Type-3 map (Planck Collaboration XIII 2014), which is defined by assuming spatially constant CO J = 2 → 1/J = 1 → 0 and J = 3 → 2/J = 1 → 0 line ratios. In the current analysis, we only assume that the same approximation holds between detectors within single frequencies, i.e., that the ratio between the CO signals observed by the 100-ds1 and 100-ds2 detectors is spatially constant. These assumptions are not strictly correct in either case, both because of real variations in local physical properties such as opacity and temperature, and because of non-zero velocities of molecular clouds that either red- or blueshift the intrinsic line frequency. Since the derivative of the bandpass profile evaluated at the line frequency also varies between detectors, the effective observed line ratio along a given line of sight also varies on the sky. As we shall see in Sect. 5, this effect represents the dominant residual systematic in a few of our frequency channels after component separation.

Thermal Sunyaev-Zeldovich – The last of the main astrophysical components included in the following temperature analysis is the thermal Sunyaev-Zeldovich (SZ) effect, which is caused by CMB photons scattering off hot electrons in clusters (Sunyaev & Zeldovich 1972). After such scattering, the effective spectrum no longer follows a perfect blackbody, but is rather given by the expression⁹ listed in Table 4. The only free parameter for this effect is the Compton parameter, y_sz, which for our purposes acts as a simple amplitude parameter. We note that the effective SZ spectrum is negative below 217 GHz and positive above this frequency, and this feature provides a unique observational signature. Still, the effect is weak for all but the very brightest clusters on the sky, and the y_sz map is therefore particularly sensitive to both modelling and systematic errors. In this paper, we only fit for the thermal SZ effect in two separate regions around the Coma and Virgo clusters, which are by far the two strongest SZ objects on the sky, in order to prevent these from contaminating the other components. While the SZ decrement for the Coma cluster is as large as −400 μK below 100 GHz on an angular scale of a few arcminutes (Battistelli et al. 2003; Planck Collaboration Int. X 2013), and for the Virgo cluster a few tens of microkelvin on a few degrees scale (Diego & Ascasibar 2008), it is only a few microkelvin for most other objects after convolution with the large beams considered in this paper. Further, full-sky SZ reconstruction within the present global analysis framework requires significantly better control of systematic effects than what is achieved in the current analysis, in particular at high frequencies. For dedicated full-sky SZ reconstruction, see Planck Collaboration XXI (2014) and Planck Collaboration XXII (2016).

Monopoles and dipoles – In addition to the above astrophysical components, the microwave sky exhibits important signal contributions in the form of monopoles and dipoles. The prime monopole examples are the CMB monopole of 2.7255 K and the average CIB amplitude (Planck Collaboration XXX 2014). The main dipole contribution also comes from the CMB, which has an amplitude of 3.3655 (3.3640) K as measured by LFI (HFI). The difference between the LFI and HFI measurements of 1.5μK is within quoted uncertainties (Planck Collaboration I 2016).

Ideally, the CMB dipole contribution should be removed during map making (Planck Collaboration VI 2016; Planck Collaboration VIII 2016), but because the estimated dipole has a non-zero uncertainty, and different dipoles are subtracted from the Planck and WMAP sky maps, it is necessary to account for residual dipoles in each individual map. In this paper, we adopt as delta function priors the nominal Planck dipole estimates for a small subset of reference frequencies (namely the 30, 100-ds1, 143-ds1, 545, and 857 channels), as well as the parameters derived by Wehus et al. (2016) for the 408 MHz map. That is, the dipoles of these six channels are fixed at their nominal values, and, together, they anchor the dipole solutions for all astrophysical component maps.

Likewise, we adopt the nominal monopole parameters for a subset of frequency channels, anchoring the effective offsets of each component to external values. Specifically, we adopt the HFI CIB monopole predictions, as listed in Table 5 of Planck Collaboration VIII (2016), for the 100-ds1, 143-ds1 and 545-2 monopoles, and the zero-points derived by Wehus et al. (2016) for the WMAP Ka and 408 MHz bands.

Thus the following analysis does not derive independent absolute estimates for either monopoles or dipoles, but relies critically on external priors for these values. If, say, the CIB monopole prediction is significantly updated through additional observations, this error will translate directly into an error in the zero level of our thermal dust map. On the other hand, our analysis does provide an independent internal consistency check among the adopted priors, and, as we shall see, no major anomalies are found.

Other components – Finally, we mention briefly some other sub-dominant effects that are either neglected or only indirectly accounted for in our model. First, similar to the thermal SZ effect, the kinetic Sunyaev-Zeldovich effect is caused by CMB photons scattering on hot electrons (e.g., Rephaeli 1995; Planck Collaboration XXII 2016, and references therein). However, in this case a non-zero bulk velocity in the electron population is the defining feature. When interacting with moving electrons, the photons effectively receive a Doppler shift proportional to the bulk velocity. In this process, the CMB blackbody spectrum is conserved, but its temperature is slightly changed. The kSZ effect is therefore fully degenerate with the primordial CMB anisotropies when considered pixel-by-pixel, and can only be disentangled using spatial information and cross-correlations. However, this effect is small, and we neglect it in the rest of the paper.

Second, the cosmic infrared background Planck Collaboration XVIII (2011; Planck Collaboration XXX 2014) consists of redshifted thermal dust radiation from distant galaxies. By virtue of being a cosmological signal, it is statistically isotropic in the sky, and can therefore be described analogously to the CMB in terms of a sum of a constant offset at each frequency plus small fluctuations imprinted on top, tracing the column density of the emitting matter integrated over redshift. In this paper, we account for the mean temperature through the monopole parameter of each frequency map, but neglect the spatially varying component. Since the CIB frequency spectrum is very similar to the Galactic thermal dust spectrum, with typical mean greybody parameters of β_CIB = 1.4 ± 0.2 and T_CIB = (13.6 ± 1.5) K (Gispert et al. 2000), these fluctuations are effectively absorbed by our thermal dust model. Future work will attempt to break this degeneracy by exploiting the different spatial power spectra of the two components, and through the use of external priors.

Third, there has been some discussion in the literature concerning the existence of magnetic dipole emission from dust grains (Draine & Lazarian 1999; Draine & Hensley 2013), with tentative evidence reported in Planck Collaboration Int. XXII (2015). Suggestions have even been made that such emission might contaminate previously published CMB maps at significant levels (Liu et al. 2014). Planck Collaboration XXV (2016) comments briefly on this question, while in this paper we simply observe that the current model is able to reproduce the data with statistically acceptable accuracy without such a component, and it is neglected in the following baseline analysis.

Fourth and finally, extra-Galactic point sources contribute significant power at both low and high frequencies (Planck Collaboration XXV 2016). In the following analysis, in which we fit for foreground parameters pixel-by-pixel, these sources end up in the component map with the most similar frequency spectrum, for instance the free-free or thermal dust emission maps. In an attempt to minimize such contamination, we preprocessed each input map by subtracting catalogues of known resolved sources (Bennett et al. 2013; Planck Collaboration XXVIII 2014) before performing the analysis. However, due to source variability, beam asymmetries, and uncertainties in the actual catalogues, we invariably found worse results with such pre-processing. A more promising approach is to fit for each source jointly with the diffuse components, as for instance described by Bull et al. (2015). This, however, is left for future work.

4.2. Sky components in polarization

Whereas the microwave sky in total intensity exhibits very rich astrophysics, in principle requiring between 10 and 20 different physical components for a proper model (depending on the level of detail required), the same is not true for polarization. At the sensitivity level of current experiments, only three diffuse components have been clearly detected, namely CMB, synchrotron, and thermal dust.

CMB – Thomson scattering between electrons and photons is an intrinsically polarized process. If the incident radiation surrounding a given electron is fully isotropic, no net polarization will be observed from that electron. However, if the radiation forms a local quadrupole, a net non-zero CMB polarization signal results (Zaldarriaga & Seljak 1997; Kamionkowski et al. 1997; Hu & White 1997). Thus, any process that generates quadrupolar structures at the time of recombination or during the epoch of reionization will result in a CMB polarization signal. The scalar fluctuations that are responsible for creating the CMB temperature anisotropies will form so-called “E-mode” patterns, in which the polarization direction is either orthogonal or parallel to the wave direction. Inflationary gravitational waves and weak gravitational lensing of CMB E-modes produce so-called “B-mode” patterns, in which the polarization direction is rotated −45deg or + 45deg with respect to the wave direction. These two are the only known sources of such B-modes of primordial or high redshift origin.

Synchrotron – Radiation emitted from relativistic electrons spiraling in a magnetic field is intrinsically highly polarized, with about seven times more energy being emitted in the plane of the electron’s motion than in the orthogonal direction. For a perfectly regular magnetic field, the synchrotron polarization fraction may exceed 70% (Pacholczyk 1970), although for realistic fields it is usually significantly less. WMAP found a polarization fraction of about 3% in the Galactic plane and about 20% at high Galactic latitudes (Page et al. 2007). Other analyses report polarization fractions as high as 50% on large angular scales (Kogut et al. 2007; Vidal et al. 2015). The new Planck measurements at 30 and 44 GHz complement these observations, and allow us to put independent constraints on this quantity (Planck Collaboration XXV 2016).

Thermal dust – Magnetic fields also have an important effect on aspherical dust grains, in the sense that they tend to align the grains’ major axes with the local magnetic field direction. This alignment naturally results in net microwave polarization with a thermal dust spectrum. Furthermore, with the frequency coverage and sensitivity of the HFI channels, Planckis ideally suited to measure this signal with very high accuracy, as has already been demonstrated through a series of recent papers, including Planck Collaboration Int. XIX (2015), Planck Collaboration Int. XX (2015), Planck Collaboration Int. XXI (2015), Planck Collaboration Int. XXII (2015), and Planck Collaboration Int. XXX (2016). One important goal of the current paper is to provide direct access to these observations in a computationally convenient form.

Neglected components – Considering the other components that are relevant for intensity analysis, we first note that the free-free emission process is intrinsically independent of direction, and therefore naturally unpolarized, although some polarization may arise due to Thomson scattering on electrons near the edges of strong H ii regions. This effect, however, is smaller than 10% for edges of optically thick clouds (Keating et al. 1998), and negligible away from any sharp edges. Averaged over the full sky, free-free emission is observationally constrained to be less than 1% polarized (Macellari et al. 2011).

Next, spinning dust emission is generated by small dust grains, and these generally align weakly with the local magnetic field. However, rotational energy level splitting dissipates energy, and this leads to a low level of grain alignment (Lazarian & Draine 2000). As a result, the polarization fraction may be up to 1–3% between 10 and 30 GHz, but falling to less than 1% at higher frequencies. For a recent review of current constraints, see Rubiño-Martín et al. (2012).

Third, polarized CO emission was first detected by Greaves et al. (1999) at the roughly 1% level near the Galactic centre and in molecular clouds. This signal may therefore in principle be detected in the very sensitive Planck 100 and 217 GHz frequency channels. However, as described in Planck Collaboration VIII (2016), one of the most important systematic effects in the Planck polarization data set is temperature-to-polarization leakage, one component of which is precisely leakage of the CO signal. For now, therefore, we conservatively mask any regions with strong CO intensity detections from all polarization maps (see Sect. 6 for further discussion). This issue will be revisited in the future.

Contributions from the CIB, extragalactic point sources (Tucci & Toffolatti 2012), the Sunyaev-Zeldovich effect, etc., are all small compared to the noise level of Planck, and neglected in the following.

4.3. Instrumental effects

In addition to astrophysical components, our model includes two important instrumental effects, relative calibration between detectors and bandpass uncertainties.

Relative calibration – A calibration error scales a map up or down in amplitude, affecting all components equally. For astrophysical foregrounds, such errors are therefore strongly degenerate with spectral indices. The CMB, however, is a nearly perfect blackbody across the relevant frequency range, and provides a powerful relative calibration source. In the following, we (arbitrarily) adopt the 143-ds1 channel as our absolute calibration reference. For any channel between 30 and 353 GHz, we then use CMB fluctuations between multipoles ℓ_min = 25 and 100 to determine the relative calibration by means of a simple cross-correlation, (15)where s_cmb = ∑ _ℓ,ma_ℓmY_ℓm is the usual spherical harmonic transform. This equation effectively replaces Eq. (10) in the Gibbs sampling chain for the relevant channels. The cross-correlation coefficient is evaluated over pixels admitted by the PM61 processing mask (Fig. 3), and the multipole range 25 ≤ ℓ ≤ 100 is choosen to minimize contamination from diffuse foregrounds on the low end and confusion with beam uncertainties on the high end.

Strictly speaking, this estimation step is a violation of the Gibbs sampling algorithm as defined in Eq. (10), since that formally should include all signal components over the full sky, not just the CMB component evaluated over parts of the sky. However, in this special case we trade mathematical rigour for the sake of increased robustness with respect to systematic effects.

The 545 and 857 GHz channels are treated as special cases. As discussed in Sect. 3, these channels are significantly affected by far sidelobe and destriping error contamination. We therefore exclude the 545-1, 857-1, 857-3, and 857-4 channels entirely, and add white regularization noise of 0.01 MJy sr^-1 to the remaining channels. Unfortunately, this noise addition has the side-effect that any CMB fluctuations at 545 GHz are obscured by random white noise, to the extent that it is no longer possible to calibrate the 545 GHz channel with CMB fluctuations within the usual Gibbs chain. As a partial solution to this problem, we perform a dedicated pre-analysis that is identical to the main analysis, except that no regularization noise is added to the high-frequency channels. The 545-2 calibration is then fixed at the resulting value in the subsequent main analysis, while the other 545-4 and 857-2 calibrations are fitted using the full Gibbs expression, including full-sky (up to the PM99.6 processing mask) and all-component observations. We estimate that the uncertainty on the 545-1 recalibration factor derived in this paper is 1–2%, which translates into a 3–6% uncertainty on the 857 GHz calibration through the thermal dust scaling. Note, however, that the 857 GHz calibration is almost perfectly correlated with the 545 GHz calibration. These two values represent the most significant sources of systematic error in the entire temperature analysis.

Finally, the calibration of the two lowest frequencies (408 MHz and WMAP K-band) are fixed to their nominal values. The estimated calibration uncertainty of the 408 MHz channel is 10% (Haslam et al. 1982), and as a result our synchrotron model also has a 10% calibration uncertainty. We note that the 408 MHz calibration is almost perfectly degenerate with the GALPROP scale frequency, , a degeneracy that is only partially broken by the presence of free-free emission. For the WMAP K-band, the main problem lies in the considerable model uncertainty of the spinning dust component, which prevents an effective CMB cross-correlation at the required sub-percent accuracy level. For now, we consider the default CMB dipole-based K-band calibration to be more robust than a CMB cross-correlation-based estimate.

Bandpass errors – The bandpass profile of a detector provides the information required to convert between a given foreground spectrum and the actually measured signal, both in terms of unit conversion and effective colour corrections (Planck Collaboration V 2014; Planck Collaboration IX 2014). Therefore, a measurement error in the bandpass profile essentially translates into a spectral-index-dependent multiplicative scale error. As such, bandpass measurement errors are almost perfectly degenerate with calibration errors, except in one crucial aspect. For maps that are calibrated on the CMB dipole, the amplitude of the CMB fluctuations is fully independent of any bandpass assumptions. This fact holds the key to separating bandpass from calibration effects. While calibration errors affect both foreground and CMB components equally, bandpass measurement errors only affect foreground components.

As shown in the following, bandpass measurement uncertainties are important for both LFI and HFI in high-precision component separation. Figure 5 provides a simple illustration of this, in the form of the fractional difference between the 143-ds1 and 143-ds2 detector maps, (16)The top panel shows the fractional difference map computed from the observed Planck 2015 data and smoothed to 40′ FWHM, and the bottom panel shows the same from the 2015 FFP8 simulation (Planck Collaboration XII 2016), which is based on the nominal bandpass profiles. The effective frequencies for these two channels, based on laboratory measurements, are listed in Table 4 of Planck Collaboration IX (2014), and are (144.22 ± 0.02) and (145.05 ± 0.02) GHz for a spectral index of α = 2. The nominal 143-ds2 effective frequency is therefore 0.8 GHz higher than the corresponding 143-ds1 effective frequency, and this holds almost independently of the assumed spectral index. Given these numbers, and recognizing that the 143 GHz channel is dominated by thermal dust, which has an effective spectral index of β ≈ 1.4 at this frequency, one would expect the above fractional ratio to be (1−(145.05 / 144.22)^1.4) ≈ −0.8%, as is indeed observed in the simulation. However, in the real data, we see that the fractional ratio is approximately +0.2%, strongly suggesting that 143-ds1 in fact has a higher effective frequency than 143-ds2, not lower.

An alternative explanation for this discrepancy between the model and the data might be relative calibration errors. However, a 1% relative calibration error in the 143 GHz channel is strongly ruled out by CMB cross-correlations, since the internal relative calibration of this channel is accurate to about 0.1% (see Sect. 5, Planck Collaboration I 2016, and Planck Collaboration VIII 2016). A second possibility is that the 143 GHz channel might be dominated by some other component than thermal dust, such as free-free emission. However, such a hypothesis can only explain the observed difference for one particular detector combination, while similar differences are observed between many detector maps. Furthermore, it does not match the morphology of the observed fractional difference map. The only viable explanation we have found is indeed bandpass errors with magnitudes at the ≲1 GHz level. It is worth noting that significant discrepancies among the 143 GHz detector bandpasses were noted already in Planck Collaboration IX (2014), and these qualitative findings are therefore not new.

Fortunately, there is sufficient information in the current data to self-consistently constrain and mitigate these errors. First, by employing detector-set maps instead of full-frequency maps, bandpass errors within frequency bands can be constrained with high signal-to-noise levels, as illustrated in Fig. 5. Second, because the true sky contains multiple foreground components with varying spectral indices, while there is only one true bandpass per detector that applies to all components and all pixels on the sky, even the mean bandpass correction per frequency is constrained to some level, although this dependency is conditional on the overall foreground model. Therefore, the bandpass corrections presented in the following are strongly constrained with respect to relative bandpass shifts within frequency bands, but uncertain with respect to frequency shifts across frequency bands. Intuitively, it is easy to conclude with high significance that the relative difference between the 143-ds1 and 143-ds2 effective frequencies is +0.2 GHz, not −0.8 GHz, but it is difficult to say with confidence whether both should be shifted jointly by an additional 0.2 or −0.2 GHz. That depends on the spectral index of thermal dust emission, and consequently on the 545 and 857 GHz calibration factors, which are still associated with significant uncertainties.

As for an explicit parameterization of the bandpass errors, we adopt in this paper a simple and purely phenomenological model given by a rigid bandpass translation in frequency space: (17)Here b₀(ν) denotes the nominal bandpass, b(ν) is the fitted bandpass, and Δ_ν is the effective frequency shift between the two. We emphasize that this is not to be taken as a physical model of the true bandpass errors, and more realistic models should include information specifically regarding the tails of the bandpasses, the separation between main bandpass modes, and overall large-scale tilts, not only centre frequency shifts.

5. Temperature analysis

5.1. Model motivation

Our baseline Planck 2015 astrophysical model, derived from a joint analysis of the Planck, WMAP, and 408 MHz observations, includes the following components: (1) CMB; (2) synchrotron emission; (3) free-free emission; (4) two spinning dust emission components; (5) CO line emission at 100, 217, and 353 GHz; (6) 94/100 GHz line emission; (7) thermal dust emission; and (8) thermal SZ emission around the Coma and Virgo clusters. In addition, the model includes monopole, dipole, recalibration, and bandpass corrections for each channel.

Many model variations were considered in the preparatory stages of this analysis, including: (1) broken or smoothly-curved synchrotron frequency spectra; (2) free-free models with fixed power-law indices; (3) one-component spinning dust models, based both on physical SpDust templates and phenomenological second-order polynomial frequency spectra; (4) two-component thermal dust models; (5) thermal dust models with flattening at low frequencies; (6) Type-3 CO line emission; (7) no 94/100 GHz line emission component; (8) full-sky SZ reconstruction; (9) no bandpass corrections; and (10) various priors for relevant components. It is impossible to fully describe all of these topics within a single paper, and we do not attempt to do so. Some low-frequency model variations are, however, discussed in Planck Collaboration XXV (2016).

In this paper our focus is the baseline model, which is the simplest sufficient model considered to date, as defined in terms of three criteria. First and foremost, the baseline model has to provide a statistically acceptable fit to the observations over the bulk of the sky, including most of the Galactic plane. Second, it should have the smallest possible number of free degrees of freedom that still result in a statistically acceptable fit to the data. Third and finally, in cases where there are free choices, physically well-motivated models are preferred over phenomenological models.

Three examples illustrate our approach. First, for synchrotron emission we have found that a broken power-law model, i.e., a power-law model with fixed but different spectral indices above and below some break frequency (say, 3 GHz) provides an equally good χ² fit as the GALPROP model do. A straight power-law spectrum, on the other hand, does not fit the data, because the 408 MHz map appears dimmer than expected from a straight extrapolation from 23 GHz to 408 MHz, assuming the spectral index of β ≈ −3.0 to −3.2 required to fit higher frequencies. The reasons for preferring the GALPROP model are therefore not data-driven, but rather that it requires less tuning (e.g., no choice of break frequency, no free spectral indices, very weak dependency on ), and that it is based on a physically well-motivated calculation (Strong et al. 2011; Orlando & Strong 2013). The cost, however, is less flexibility for tracing real spatial variations in the synchrotron spectral index, and thereby possibly greater cross-talk between synchrotron, free-free, and spinning dust. Nevertheless, in the absence of sufficiently strong data to disentangle these variations, we consider this a lesser evil than introducing a very strong statistical degeneracy between synchrotron, free-free, and spinning dust.

A second important example is the spinning dust component, which is currently implemented in terms of two independent SpDust components, one with a free peak frequency, ν_p, per pixel, and one with a spatially constant peak frequency, for a total of three free effective spinning dust parameters per pixel. While the introduction of the first component is unavoidable when combining the Planck and WMAP observations, the necessity of the second is more subtle. Without it we invariably find highly significant residuals (many tens to a few hundreds of microkelvin in the Galactic plane) in the 60–70 GHz frequency range, with dust-like morphology. This strongly suggests a model problem with either spinning or thermal dust (or both), but so far the only acceptable solution we have found is effectively to “widen” the SpDust spectrum slightly, which is most easily implemented by adding a second independent component. Introducing, say, a flatter thermal dust index at around 100 GHz tends to exacerbate rather than ameliorate the problem. We emphasize, however, that the current two-component spinning dust model is a purely phenomenological model introduced in the absence of physically well-motivated alternatives. We fully anticipate that additional low-frequency data or further theoretical developments will improve this model significantly in the future.

Third, we currently adopt a one-component greybody model with free emissivity index and temperature for thermal dust emission. Experimenting with various two-component alternatives, we find only one absolute requirement on the thermal dust model for frequencies up to 857 GHz, namely that at least three free thermal dust parameters per pixel are required to achieve an acceptable fit to the high-frequency observations. Whether those are { A_d,β_d,T_d }, , or even (!) is not clear from the current data set¹⁰. On the other hand, it is clear that additional parameters are not required to model the current data set. For now, we prefer the one-component model simply because it has fewer global parameters than a corresponding two-parameter model, i.e., it requires no spatially fixed β_d and T_d for a second component, and therefore requires less tuning. This is not to be taken as a statement on the relative physical merits of the two models, however. Additional high-frequency observations are required to distinguish between them.

5.2. Data preview

Before presenting the actual baseline results, it is useful to visually consider a few specific data combinations in order to gain some intuition regarding the power of these data for constraining specific parameters. In Fig. 6, we show four different sky maps, each of which highlights an important and distinct feature in the data. A similar discussion based on an internal linear combination (ILC) method is presented in Sect. 3 of Planck Collaboration XXV (2016).

Starting with the top panel, we plot the effective power-law index of the combined Planck, WMAP, and 408 MHz data, when fitting only a single power-law model at low-frequencies, as opposed to fitting individual synchrotron, free-free, and spinning dust components. All other components are defined by the usual baseline model. On the low-frequency side, this approach is thus identical to the 2013 model (Planck Collaboration XII 2014), although the data volume and frequency range are significantly increased. Despite these differences, the 2013 and 2015 low-frequency spectral index maps agree reasonably well, with a mean and standard deviation difference of 0.16 ± 0.14.

The main features in this spectral index map can be broadly characterized into two types of spatial signatures. First, we see many distinct regions (e.g., the Gum nebula, Orion, Zeta Oph, and the Cygnus complex) with a shallow spectral index of β ≈ −2 and morphology as expected from various free-free tracers (e.g., Alves et al. 2015). Second, there are extended dark regions surrounding the Galactic plane, with very steep spectral indices of β ≲ −3.6, and morphology similar to that of thermal dust emission. In addition, it is possible to see some weak hints of the North Galactic Spur, a strong synchrotron emission region, but since we adopted a synchrotron-like prior of β ~ N(−3,0.2) for the low-frequency power-law index in this analysis, it is difficult to distinguish this region from the prior-dominated background. The main point, however, is that even with such minimal modelling, there is strong evidence for at least three distinct physical components at low frequencies, namely: synchrotron, with β ≈ −3; free-free, with β ≈ −2; and a dust-like component, with β ≲ −3.6. Among all the spectra shown in Fig. 4, only a spinning dust component with ν_p ≲ 30 GHz can reasonably account for the latter.

In the second panel we show the ratio between the 857 and 545 GHz frequency maps, masking all pixels for which the 545 GHz amplitude is smaller than ten times its monopole. The spatial variations seen in this map cannot be explained either in terms of calibration or bandpass errors (because it is a ratio map, and either of those errors primarily affects the scale, not the morphology) or in terms of absolute offsets (because of the high mask threshold). They are robust features of the Galactic thermal dust emission, and must be explained by any Galactic model that include these observations. In fact, within our baseline model, these features can only be explained in terms of a spatially varying dust temperature. To be specific, the thermal dust appears hotter (i.e., has larger 857-to-545 ratio; see Fig. 4) near the Galactic centre than in the outer Galaxy (quadrants 2 and 3). Adopting a spatially constant dust temperature is no longer possible, and the only reason that the corresponding 2013 analysis could produce meaningful results with such an assumption was that it considered frequencies only up to 353 GHz, and also focused primarily on high Galactic latitudes (Planck Collaboration XII 2014).

The third panel shows the straight difference between the 100-ds1 and 100-ds2 maps. Because of the very short lever arm between these two frequencies, essentially all continuous emission mechanisms cancel in this difference, leaving only the sharp line emission from CO J = 1 → 0, as well as instrumental noise. This demonstrates the power of employing detector maps in the following analysis as opposed to frequency maps; using a fine-grained data set vastly increases our ability to extract line emission.

Whereas CO line emission was studied extensively in the Planck 2013 release, an additional new 94–100 GHz line feature is included in the current release (see Sect. 4.1). The primary contribution to this component is visualized in the bottom panel of Fig. 6, in terms of the WMAP W3−W2 difference map. The small but bright features near the Galactic centre (the maximum amplitude is 740μK) are line emission within the W-band, as discussed in Sect. 4.1. This line is also covered by the two Planck 100 GHz channels. To prevent this additional emission from contaminating the main CO estimates, and to achieve an acceptable fit for the WMAP W-channels, we fit for this additional line emission component at 94/100 GHz.

5.3. Baseline model

We are now ready to present the Planck 2015 baseline diffuse astrophysical component model in temperature, as derived from the Planck, WMAP, and 408 MHz observations. An overview of the data products delivered through the Planck Legacy Archive (PLA)¹¹ is given in Table 5, including file names and FITS extension numbering, as well as summary statistics in the form of posterior maximum, mean, and rms values.

We start our review by considering data that are smoothed to a common resolution of 1° FWHM, representing the most complete data set and model possible, given the current data set. Processing these observations through the analysis pipeline outlined in Sect. 2, and adopting the baseline model described in Sect. 4.1, we obtain the set of maximum-posterior astrophysical parameter maps shown in the top panels of Figs. 7–19. The bottom panels show the corresponding rms maps found by computing the standard deviation over the ensemble of Gibbs samples. Please note that these rms maps account only for statistical errors from instrumental noise, and not for systematic uncertainties due to modelling errors. They are therefore only meaningful for pixels for which the goodness-of-fit is acceptable. Most importantly, they do not represent the true errors in the Galactic plane, where modelling errors dominate statistical errors.

Instrumental parameters, as well as monopole and dipole coefficients, are listed in Table 6, and Fig. 20 provides a visual comparison of the calibration factors for each Planck and WMAP detector map¹². Global (i.e., spatially uniform) astrophysical parameters are listed in Table 7. The full marginal uncertainties of these parameters are dominated by modelling, not statistical errors. For this reason, the tabulated uncertainties are computed through realistic end-to-end simulations, as implemented in the Planck 2015 FFP8 simulations (Planck Collaboration XII 2016). These simulations were analysed blindly as part of the CMB validation efforts (Planck Collaboration IX 2016), using the exact same machinery as used in this paper. Further, they are based on a significantly different foreground model than the baseline model adopted here, and they therefore at least partially account for modelling errors, as well as known systematic and mapmaking uncertainties. The only differences in the fitted model compared to the present baseline are that it includes only one spinning dust component and no 94/100 GHz line emission or SZ components.

Based on these simulations, we estimate the uncertainties on the calibration and bandpass measurements directly from the simulations, comparing output parameters against known inputs. The monopole and dipole uncertainties are, however, first based on the statistical errors derived from the Gibbs chains, and then rounded up (where necessary) to be no smaller than the corresponding FFP8 simulation uncertainties¹³. Thus they correspond to the maximum of the instrumental and the modelling errors. Furthermore, we reemphasize that the reported monopole and dipole values are conditional with respect to the nominal Planck, H i and synchrotron priors, as defined by Planck Collaboration VI (2016), Planck Collaboration VIII (2016), and Wehus et al. (2016).

The uncertainties in the CO line ratios listed in Table 7, which are model-dominated, are also taken from the FFP8 simulations. The 94/100 GHz line ratios, however, are noise dominated, and they are therefore taken directly as the standard deviation of the posterior distribution, i.e., from the Gibbs samples. Finally, we do not quote any uncertainties on the peak frequency of the second spinning dust component or the frequency scale factor of the synchrotron GALPROP model. Both of these are completely dependent on other parameters in the model, and the full joint distribution is highly non-Gaussian. As a result, we only quote the maximum-likelihood values for these, and notice that their full marginal uncertainties are very large, possibly up to 50% or more.

5.3.1. Goodness-of-fit

Next, the goodness-of-fit of the model is summarized in Figs. 21–22 and in Table 8. Starting with Fig. 21, residual maps for each channel are shown on the form (d_ν−s_ν), and these give the most complete view of the performance of the model out of all the statistics considered in the following. In addition, they provide a direct visual summary of remaining low-level systematics in the data, which should prove useful for future reprocessing of the same data. Overall, we see in these maps that the model provides an excellent fit to all channels at high Galactic latitudes. Please note that the colour scales vary from ± 2 μK at the CMB dominated Planck HFI channels, through ± 10 μK for most of the WMAP channels, to ± 1K and ± 0.05 MJy sr^-1 for the 408 MHz and 545 GHz channels. In fact, the solution is sufficiently free of artifacts that the dominant systematic effect in the 100-ds2 residual map is Galactic rotation projected into the CMB frequencies by the CO J = 1 → 0 emission line. The same feature is also seen in several other channels, although at slightly lower significance.

The top panel of Fig. 22 shows the corresponding χ² map as defined by Eq. (14), but summed only over the accepted frequency channels. The bottom panel shows a histogram of these χ² values including all pixels admitted by the PM61 processing mask.

Together, these figures provide a useful qualitative summary of the goodness-of-fit of the baseline model. However, providing a corresponding rigorous statistical description is complicated, because the effective number of degrees of freedom per pixel is not well defined. The usual approach of subtracting the number of free parameters from the number of data points is not applicable, because informative priors (most notably the positivity prior) eliminate large parts of the parameter space. A parameter that is prior-dominated therefore does not reduce the number of degrees of freedom by unity, but only a fraction of unity. Second, smoothing the data to a common resolution of 1deg FWHM introduces noise correlations between pixels, and these are not accounted for in the noise description. The overall χ² distribution will therefore be broader than a corresponding distribution with no correlated noise. To estimate the effective number of degrees of freedom, we therefore fit a scaled χ² distribution to the empirical χ² distribution, including only the very cleanest parts of the sky, where actual foreground contamination is minimal. We adopt the conservative PM61 mask for this. The resulting best-fit χ² distribution reads (18)and we accordingly estimate that the empirical number of degrees of freedom is 20, and the correlated noise scaling factor is 1.37. The former of these suggests that our model effectively contains 32−20 = 12 free parameters, not 14 as obtained by naively counting free parameters per pixel. In other words, the combined effect of all priors is to remove 2 of 14 degrees of freedom, indicating that the model is indeed highly data driven. he latter number suggests that the white noise approximation underestimates the true noise level by 37%. The 99% confidence χ² range from the analytic fit is 11 <χ²< 59.

The middle panel of Fig. 22 shows the mask obtained by thresholding the χ² map smoothed with a 1° Gaussian kernel at a value of 50.¹⁴ This mask is called the 93% likelihood mask (LM93), and constitutes the primary confidence mask for the Planck 2015 model. Also, together with the CMB solution in Fig. 7, this mask defines the main inputs to the Planck 2015 low-ℓ CMB temperature likelihood (Planck Collaboration XI 2016). We note that this mask removes many bright extra-galactic point sources at high Galactic latitudes; the algorithm adopted in this paper treats point sources and diffuse emission identically through pixel-by-pixel fits, and any subsequent analysis of the resulting component maps should take these sources into account either through explicit masking, as done here, or by post-processing fits.

The second column of Table 8 lists the rms of each residual map outside the LM93 mask. The third column lists the ratio between these rms values and the corresponding instrumental noise rms, as listed in Table 1. The fourth column shows the monopole and dipole corrected median fractional residual, (d_ν−s_ν) / (d_ν−T_νm_ν), this time evaluated inside the LM93 mask. As seen in the third column, the rms values of the residual map are close to the instrumental noise for most channels, indicating an excellent goodness-of-fit, not only in terms the integrated χ², but also channel-by-channel.

In accordance with the above discussion, we see that a number of channels have residuals that are lower than the instrumental noise, with the 408 MHz and 857-2 channels being the two most striking examples. This is expected, because the χ² only has 20 degrees of freedom, whereas there are 32 individual frequency channels; the normalized mean residuals must therefore sum to less than one per channel. However, a ratio much smaller than unity typically indicates that the corresponding channel has a much higher effective signal-to-noise ratio with respect to some signal parameter than all other channels combined. In our case, the 408 MHz and 857-2 channels strongly dominate the synchrotron and thermal dust amplitudes, respectively. At first sight, one might suspect these values to indicate the presence of worrisome parameter degeneracies, which typically also can result in residuals lower than the instrumental noise. However, from the parameter maps shown in Figs. 7–19, it is visually obvious that the synchrotron and thermal dust emission amplitude maps are not degenerate with any components. Rather than degeneracies, these low residual values indicate that the current data set is non-redundant with respect to these two amplitude maps; the 408 MHz map determines almost exclusively the synchrotron amplitude, and the 857-2 map determines almost exclusively the thermal dust amplitude. The main problem with these low residuals is therefore not degeneracies, but rather lack of robustness with respect to systematics; any systematic error that may be present in the 408 MHz and 857-2 channels can and will propagate into the respective foreground amplitude maps. In order to improve on this situation in the future, recovering the currently systematics contaminated 857-1, 857-3, and 857-4 channels is imperative on the high-frequency side, and incorporating additional low-frequency observations (between, say, 1 and 20 GHz) is critical on the low-frequency side.

A similar effect is seen for a number of other channels, if not equally strongly. For instance, we see that the WMAP K-band and Planck 30 GHz channels have rms ratios of 0.38 and 0.55, respectively, and these dominate the two spinning dust amplitudes, A_sd1 and A_sd2. The 100-ds2 channel has a ratio of 0.76, and defines the CO J = 1 → 0 amplitude together with 100-ds1. Finally, the 353-1 channel has a ratio of 0.60, and this channel has the greatest pull on the dust emissivity index, β_d.

The single most important conclusion from Table 8, however, is that the baseline Planck temperature sky model is an excellent fit to the observed data, in agreement with the visual impression of Fig. 21. The residuals are largely consistent with instrumental noise over 93% of the sky, and the median fractional residual in the complementary 7% of the sky is; smaller than 0.2% for all HFI channels; smaller than 1% for all LFI channels; and smaller than 2% for all WMAP channels.

5.4. High-resolution component maps

We now consider the high-resolution intensity maps derived using the same pipeline as above, but with a reduced data set and astrophysical model. Specifically, we only include channels from 143 GHz and above, all smoothed to a common resolution of FWHM. The signal model includes CMB, thermal dust, and CO J = 2 → 1 and 3 → 2 emission lines coadded into one map, similar to the Type-3 map in the 2013 data release¹⁵. We fix all global instrumental parameters on the values listed in Table 6, and the thermal dust temperature, T_d, to the low-resolution solution, upgraded in harmonic space (to avoid pixelization effects) to a HEALPix resolution of N_side = 2048. The only free spectral parameter per pixel is now the thermal dust emissivity index, β_d.

The resulting amplitude maps are shown in the top panels of Figs. 23–25, while the bottom panels show the so-called half-ring half-difference maps, i.e., half-difference between two maps derived from independent half-ring maps (Planck Collaboration VI 2016; Planck Collaboration VIII 2016); these provide a direct estimate of the instrumental noise present in the full-mission maps. Figure 26 shows the same for the high-resolution dust spectral index.

We note that while the Galactic centre appears negative in the low-resolution CMB solution shown in Fig. 7, it appears positive in the corresponding high-resolution CMB shown in Fig. 23. This qualitative difference demonstrates the importance of modelling errors in the Galactic plane. At high resolution, our model includes only CMB, CO and thermal dust, but no dedicated low-frequency component. Any residual free-free contributions to frequencies at or above 143 GHz is therefore necessarily interpreted as a combination of CMB and CO, both of which have redder spectral shapes than thermal dust. In the low-resolution solution, on the other hand, the main problem lies in the interplay between CO and thermal dust modelling errors and high-frequency calibration uncertainties.

5.5. Comparison with independent data products

To further validate the baseline model presented in Sect. 5.3, we now compare the derived products with similar maps produced either by independent observations or through different analysis techniques, focusing in particular on spinning and thermal dust and CO emission. Synchrotron and free-free (and spinning dust) are addressed separately in a dedicated companion paper, and we refer the interested reader to Planck Collaboration XXV (2016) for full details. A short summary of that analysis includes the following main points.

• 1.

  The synchrotron estimates derived by the WMAP team (Bennett et al. 2013) using either Markov chain Monte Carlo (MCMC) or maximum entropy methods (MEM) typically have 50–70% higher amplitudes at high Galactic latitudes compared to that derived in this paper, and this increases to factors of several when including the Galactic plane. This is compensated by about twice as much spinning dust in our model compared to the WMAP models.

• 2.

  The free-free model derived in the current analysis correlates well with Hα observations. For instance, the scaling factor between the two maps in the Gum Nebula is (8.2 ± 1.3) μK R^-1. For comparison, Dickinson et al. (2003) found values ranging between 8.2 and 13.1 μK R^-1, depending on the exact position within the Gum Nebula.

• 3.

  The free-free map also shows good morphological agreement with respect to radio recombination line (RRL) observations (Alves et al. 2015), although the derived amplitude is notably higher in our map. The six brightest objects have a mean relative amplitude ratio of 1.14 ± 0.04, whereas the ten next have an amplitude ratio of 1.36 ± 0.08. Considering that RRLs are in principle a very clean probe of free-free emission, the most likely explanation for this discrepancy is leakage between synchrotron, spinning dust, CO, and free-free in the current solution. On the other hand, the derived RRL amplitudes also depends sensitively on the assumed electron temperature, and raising /Te by ≈20% would resolve much of the difference. In addition, it is worth noting that other component separation techniques, including FastMEM, CCA, and both the 9-yr WMAP MCMC and MEM analyses, all derive free-free amplitudes consistent with the Commander result (Planck Collaboration Int. XXIII 2015).

5.5.1. Dust template amplitude consistency by H I cross-correlation

We next perform an internal consistency test of the Commander dust model in the range from 23 to 353 GHz, as defined by the sum of the SpDust2 components described in Sect. 4.1 and thermal dust, by cross-correlating our dust model against GASS H i observations (McClure-Griffiths et al. 2009; Kalberla et al. 2010) covering 18% of the high Galactic latitude sky near the South Galactic pole, following the procedure of Planck Collaboration Int. XVII (2014). In particular, we compare the resulting template amplitudes against the corresponding amplitudes derived directly from cross-correlation with the raw Planck and WMAP sky maps.

Figure 27 shows the fractional difference between the resulting template amplitudes for each frequency band. Overall, the agreement is satisfactory with typically 20% differences in the 20–70 GHz range, in which spinning dust provides a larger contribution to the frequency spectrum than thermal dust. At higher frequencies, where thermal dust emission starts to dominate, the agreement improves further to around 5%, and to within 1 σ in terms of statistical uncertainties.

5.5.2. Dust SED consistency by H I and internal Planck cross-correlations

Next, we consider the robustness and consistency of the thermal dust SED model, as parametrized in terms of the two greybody parameters, β_d and T_d. Specifically, we compare the new SED estimates derived in this paper with corresponding estimates derived by H i cross-correlation at high latitudes in Planck Collaboration Int. XVII (2014), and by internal Planck cross-correlations at intermediate Galactic latitudes in Planck Collaboration Int. XXII (2015), both of which have been updated with the latest Planck 2015 sky maps.

Figure 28 compares the mean thermal dust SED derived from H i–CMB cross-correlation and the Commander estimates at high Galactic latitudes in terms of the fractional difference, (s_{H i}−s_Comm) /s_Comm. The two sets of best-fit thermal dust spectral parameters are (β_d,T_K)_Comm = (1.54,22.8 K) and (β_d,T_K)_{H i} = (1.54,21.4 K), respectively, and the two models agree point-by-point to 5–10% between 100 and 857 GHz. At 70 GHz the difference is 50%, and this is due to different spinning dust modelling; as already shown in Fig. 27, the sum of spinning and thermal dust agree to 20% in this range between the two methods. Note also that 70 GHz is very close to the foreground minimum, and these numerically large relative differences therefore correspond to small absolute differences.

This test provides a robust estimate of residual systematic errors in the Commander thermal dust model from potential zero-level and dipole uncertainties in the high-frequency HFI channels arising from zodiacal light emission and CIB residuals, as discussed in Sect. 3. Because the H i analysis is insensitive to such errors, we take the 1–2 K difference between the two as an estimate of the systematic uncertainty on the Commander thermal dust temperature at high Galactic latitudes.

At frequencies above 857 GHz we also plot the extrapolation of the new Commander model into the COBE-DIRBE wavelengths of 240, 140, and 100μm (Hauser et al. 1998). Here we clearly see that the current model breaks down beyond the Planck frequencies, with a fractional difference of 40% between the Commander model and the DIRBE 100μm observations. Including the DIRBE channels in the fit would of course reduce these fractional residuals dramatically, but only at a very significant cost of increasing the residuals at lower frequencies between 70 and 353 GHz. The simple one-component greybody thermal dust model adopted in this paper is not able to simultaneously fit this entire frequency range, because of both intrinsic complexity of the dust particle population and because of residual systematics and calibration uncertainties in the DIRBE and high-frequency Planck data. Integration of these channels requires substantial additional work, and is beyond the scope of the current paper. For a first analysis of similar type, see Planck Collaboration Int. XXIX (2016).

Next, we turn to intermediate Galactic latitudes, for which the signal-to-noise ratio is higher, but also the astrophysical composition is richer. In this case we therefore compare our results with the outputs from the internal Planck template cross-correlation analysis of Planck Collaboration Int. XXII (2015). In short, this analysis estimates the SED parameters by cross-correlating the Planck 353 GHz channel with lower frequencies over circular patches of 10deg radius. Figure 29 compares the histogram of β_d derived using this method with the corresponding Commander estimates over the same sky region. The agreement is very good, and with averages and dispersions of and , respectively. When interpreting the widths of these distributions, it is useful to return to the thermal dust spectral index maps shown in Figs. 17 and 26. These maps are quite uniform, and, indeed, at the current level of leakage between thermal dust, CO, compact objects and residual offsets, there is little convincing evidence for true spatial variation in β_d in the results presented here. If this conjecture is true, the widths of the histograms shown in Fig. 29 are primarily expressions of analysis uncertainties in the form of instrumental noise, parameter degeneracies and systematic errors, rather than true spatial variation.

5.5.3. CO line emission

Finally, we compare the CO J = 1 → 0, 2 → 1 and 3 → 2 maps derived in Sects. 5.3 and 5.4 with independently derived maps and products. As described in Planck Collaboration XIII (2014), Planck has implemented a multi-algorithm approach to CO extraction, configuring the MILCA (Hurier et al. 2013) and Commander algorithms for dedicated CO reconstruction. In 2013, this resulted in three different types of CO maps. In short, the Type-1 CO maps are built from individual bolometer maps within single frequencies, and as such are only weakly dependent on foreground extrapolations, but this insensitivity comes at a high cost in terms of instrumental noise. The Type-2 maps are built per CO line from a small sub-set of frequencies, carefully selected to be optimal for CO extraction. Since more than one frequency is involved, a more elaborate foreground model is required, such as explicit modelling of CMB, dust and free-free, although several simplifications are imposed, such as the assumption of constant dust temperature and spectral indices. Finally, the Type-3 map corresponds to a maximum signal-to-noise CO extraction in which a complete foreground model is fitted with Commander, as described in this paper, but with only a single CO amplitude per pixel and otherwise only spatially fixed line ratios accounting for scaling between frequencies.

In the present release, the MILCA-based Type-1 and Type-2 maps have been updated with the latest data, while the Commander-Ruler-based Type-3 map from 2013 has been superseded by the high-resolution Commander-only J = 2 → 1 map presented in Sect. 5.4. In addition to these high-resolution maps, we of course also provide the low-resolution line maps discussed in Sect. 5.3. Table 9 summarizes the CO-related data products provided in the current release, including angular resolution, instrumental noise, and analysis assumptions.

We start by comparing the maps derived with MILCA and Commander, both with each other and with the CO J = 1 → 0 survey presented by Dame et al. (2001). The full-sky J = 1 → 0 and 2 → 1 maps are shown in Figs. 30 and 31, while zoom-ins of the Orion region are shown in Fig. 32. All maps are smoothed to a common resolution of 15′ in these plots, except the Dame et al. survey, which has an intrinsic resolution of about 20′. For reference, the 2013 Type-3 map is shown in the bottom panel of Fig. 30.

The first three panels of Fig. 33 show T–T correlation plots between each of the three CO line maps and the Dame et al. survey, all smoothed to a common resolution of 1deg FWHM; the fourth panel shows similar correlations between the high- and low-resolution Planck products. Please note that all axes are linear, and this figure therefore highlights the very brightest CO objects on the sky. The main points to take away from these scatter plots, and the maps in Figs. 31 and 32, are the following.

• 1.

  The Type-2 and low-resolution Commander J = 1 → 0 maps agree very well internally, and also correlate strongly with the Dame et al. survey. However, they both show an overall multiplicative scaling factor of about 1.4 relative to Dame et al. This level of amplitude difference is similar to what was observed in the 2013 release (Planck Collaboration XIII 2014), and is due to a combination of bandpass uncertainties in the Planck observations and the overall 10% calibration uncertainty in the Dame et al. survey. The Type-1 J = 1 → 0 map shows bigger differences with respect to the Dame et al. survey, both in the scatter plot and the Orion zoom-in. Possible explanations include the presence of a second significant line emission mechanism, such as ¹³CO J = 1 → 0 (at 110 GHz), or, possibly, thermal dust leakage.

• 2.

  In the CO J = 2 → 1 case, the Type-2 map shows some evidence of contamination, both in the form of significant curvature in the scatter plot (top right panel of Fig. 33), and as notable diffuse emission along the Galactic plane in the Orion region and full-sky map. The agreement between the Type-1 and Commander maps is, however, good, and both show tight correlations with the Dame et al. survey.

• 3.

  For the CO J = 3 → 2 map, the scatter with respect to Dame et al. is substantial for both Type-1 and Commander, and in particular the latter shows clear evidence of curvature in the T–T plot due to contamination from thermal dust.

• 4.

  A significant residual dipole aligned with the CMB dipole may be seen in the 2013 Type-3 map, and similar residuals are also present in the 2013 Type-1 and Type-2 maps. These dipole residuals have been greatly suppressed in the new 2015 maps, due to the new estimation of the CMB dipole directly from the Planck sky maps, as described in Planck Collaboration I (2016).

Figure 34 shows histograms of the various high-resolution maps, each map being smoothed to a common resolution of 15′ FWHM and re-pixelized at a HEALPix level of N_side = 512. As already noted, here we clearly see that the Type-1 maps have significantly larger instrumental noise (wider histograms near zero) than the corresponding Type-2 and Commander maps. However, we also see that the three maps converge well at intermediate amplitudes between, say, 3 and 200 K_RJ km s^-1, for both J = 1 → 0 and 2 → 1.

Next, in Fig. 35 we compare the Type-2 and Commander J = 2 → 1 maps with targeted CO J = 1 → 0 high-latitude observations published by Hartmann et al. (1998) and Magnani et al. (2000). The same test was performed for the Planck 2013 CO maps, and full details regarding methodology and data processing can be found in Planck Collaboration XIII (2014). The top panel shows a T–T scatter plot between each of their 133 detected objects and the corresponding objects in our maps. The correlation is strong for both Type-2 and Commander, and with comparable effective line ratios in the two cases. The bottom panel shows histograms including positions for which Hartmann/Magnani et al. did not find any significant CO detections, and, correspondingly, our maps also do not exhibit any detections at these positions; both histograms are consistent with noise.

Finally, as in the 2013 analysis (Planck Collaboration XIII 2014), we have also cross-correlated the Type-1, Type-2, and Commander J = 2 → 1 maps against the CO J = 2 → 1 AMANOGAWA-2SB survey (Yoda et al. 2010; Handa et al. 2012), and find very good morphological agreement among all maps. The best-fit slopes are 0.95, 0.82 and 0.86 for the Type-1, Type-2, and Commander maps, respectively. Thus, while the Type-1 CO J = 2 → 1 map is more noisy than the other two, it is also less affected by dust and free-free contamination, and it is a robust estimation of CO J = 2 → 1 in the Galactic plane. From this test, we also estimate that the overall re-calibration factors for the Type-2 and Commander CO J = 2 → 1 maps are 15–20%.

Combining all of these results, we make the following recommendations. First, the Commander CO J = 2 → 1 map supersedes the 2013 Type-3 map as the primary “CO detection” map, and is also our preferred CO J = 2 → 1 map. The main advantages of this map as compared to the Type-2 map are higher angular resolution ( versus 15′ FWHM), slightly lower noise (0.09 versus 0.11 K_RJ km s^-1), and a tighter correlation with respect to the Dame et al. survey (see Fig. 32). However, in specific regions where the signal is such that the noise level of the Type-1 CO J = 2 → 1 map becomes less of an issue, i.e., in the Galactic plane, the Type-1 map is a better alternative in terms of overall calibration and contamination. For a dedicated J = 1 → 0 analysis, we recommend the Type-2 map, which has higher angular resolution than the corresponding Commander map (15′ versus 1deg FWHM), with similar correlation properties with respect to Dame et al. Finally, for CO J = 3 → 2 we recommend the Type-1 map, due to a significantly stronger correlation with respect to Dame et al. than the corresponding Commander map.

6. Polarization analysis

We now turn our attention to the 2015 Planck measurements of polarized foregrounds. As described in Sect. 3, and summarized in Table 1, we employ full co-added frequency maps for this analysis, in order to maximize the signal-to-noise ratio, and we fit only CMB, synchrotron, and thermal dust emission amplitudes, fixing all calibration and spectral parameters to their temperature counter-parts (see Sect. 4.2). The main data products are summarized in the bottom two sections of Table 5.

From an algorithmic point of view, the polarization analysis is essentially identical to the temperature analysis, with the same summary statistics and goodness-of-fit statistics applying equally well to the Stokes Q and U parameters as to the temperature (i.e., Stokes I). We therefore proceed in the same manner as for the temperature case, and first present the main data products, then discuss internal goodness-of-fit tests. However, in this case we additionally discuss the EE and BB angular power spectra for synchrotron and thermal dust emission, recognizing the high importance of these quantities in current cosmology. Finally, we consider external validation.

6.1. Baseline model

We start by presenting the full-sky Stokes Q and U parameter maps for polarized CMB, synchrotron, and thermal dust in Fig. 36, and a 20deg × 20deg zoom-in of the south ecliptic pole in Fig. 37; corresponding CMB maps of the north ecliptic pole, including decompositions into E and B mode polarization, are shown in Planck Collaboration IX (2016). We note that the CMB map has been high-pass filtered with a cosine-apodized harmonic space filter, removing all structures below ℓ = 20, in order to suppress large-scale instrumental systematics. Additionally, the Galactic plane has been replaced with a constrained Gaussian realization in order to avoid ringing effects during filtering (see Planck Collaboration IX 2016 for further details).

Figures 38 and 39 show the polarization amplitude and direction, as defined by direct and naive estimators, We note that the former of these is a quadratic estimator, and therefore biased by instrumental noise, and the latter is also a nonlinear function. More sophisticated approaches are applied and described in Planck Collaboration Int. XIX (2015), which properly accounts for this bias. However, in this paper the main purpose of these quantities is visual interpretation, rather than quantitative model comparisons, and the naive estimators are then useful for providing information regarding both the underlying physical structures and the instrumental noise levels. Specifically, the noise level at a given position on the sky can be inferred from the polarization amplitude maps in Fig. 38 by comparing the intensity of any region to the deep fields near the ecliptic poles, where the Planck signal-to-noise ratio is maximum.

Considering first the CMB maps in Fig. 37, it is clear that Planck measures primordial polarized anisotropies with high signal-to-noise ratio at 10′ FWHM angular scales. Furthermore, it is also visually obvious that the CMB signal is strongly dominated by primordial E-modes, in the form of a clear +-type pattern in Stokes Q and a ×-type pattern in Stokes U; see Sect. 4.2 for further discussion.

Of course, this particular region of the sky is more deeply observed than any other region on the sky (except the north ecliptic pole), and the apparent high signal-to-noise level seen in the middle of the panel is not representative of the full sky. However, due to very sharply defined features in the Planck scanning strategy in this region, the signal seen near the bottom left corner of these plots is indeed representative, within a factor of a few.

In the same region, we also clearly see the Large Magellanic Cloud (LMC) in both synchrotron and thermal dust emission (see Planck Collaboration XXV 2016 for a detailed analysis of the LMC), but only very weakly in the CMB maps. Indeed, the primary signature of the LMC is a slight increase in variance rather than a systematic bias. This is quite reassuring in terms of CMB reconstruction fidelity, since the LMC represents one of the richest astrophysical objects on the sky.

In this respect it is worth recalling that we fix all calibration and spectral parameters (thermal dust index and temperature and synchrotron spectrum) in the polarization analysis to those derived in the temperature analysis. Several attempts have been made at re-estimating these parameter independently from the polarization observations, but we find that the resulting parameters invariably become significantly biased by the same large-scale systematics that are responsible for the remaining large-scale residuals in the CMB map (Planck Collaboration II 2016; Planck Collaboration VIII 2016; Planck Collaboration IX 2016). However, other analyses that explicitly exploit spatial correlations (e.g., template fitting) to suppress such systematics have been able to produce robust results, and yield only small differences between the temperature and polarization spectral indices. For instance, Planck Collaboration Int. XXII (2015) reports full-sky estimates of β_d = 1.51 ± 0.01 for intensity and β_d = 1.59 ± 0.02 for polarization, corresponding to a difference of only 3.6σ even after averaging over most of the sky. Thus, assuming identical temperature and polarization spectral indicies is a very good approximation at the precision level of the current data, considering the additional stability with respect to instrumental systematic errors it provides.

6.1.1. Goodness-of-fit

We now consider the statistical goodness-of-fit of this simple baseline model, following the same procedure as for the temperature analysis. First, Fig. 40 shows the residual maps d_ν−s_ν, for each of the seven Planck frequency maps, all smoothed to a common resolution of 40′ FWHM angular scale. Please note that the colour range is linear between ± 5 μK, and the same in all panels. As expected, we see that 143 GHz is the most sensitive frequency channel, followed by the 100 and 217 GHz channels. In addition to instrumental noise, these channels also exhibit large-scale patterns tracing the Planck scanning strategy at the ≲0.5 μK level. Although small in an absolute sense, it is important to recall that the expected peak-to-peak amplitude of a cosmological signal from reionization corresponding to an optical depth of, say, τ ≈ 0.07 is also about 0.5 μK (Planck Collaboration XI 2016). As a result, we do not consider the CMB polarization map presented here suitable for cosmological parameter estimation on large angular scales. Instead, the Planck 2015 low-ℓ polarization likelihood relies only on the 30, 70, and 353 GHz data, for which instrumental systematics are subdominant (Planck Collaboration VI 2016; Planck Collaboration VIII 2016; Planck Collaboration XI 2016).

Next, the top panel of Fig. 41 shows the corresponding χ² map, co-adding over both frequencies and Stokes parameters. Compared to the temperature case, it is here much easier to determine the appropriate number of degrees of freedom, since no spectral parameters are fitted to the data, and no positivity priors are imposed on the amplitudes. Specifically, there are in total 14 data points (two Stokes parameters in seven frequencies) and 6 free parameters (two Stokes parameters in three components), resulting in a net 8 degrees of freedom. The nominal 95% confidence region for this number of degrees of freedom is χ² = (2,17).

As usual, the Galactic plane is the most significant feature in the χ² map. Furthermore, when comparing this χ² map (and the individual frequency residual maps) with the various component amplitude maps derived in the temperature analysis, we find strong correlations between the χ² map and the CO emission maps. This is expected from the mapmaking analyses presented in Planck Collaboration XI (2016), and, as already noted, a general recommendation regarding these maps is to reject any pixels with significant CO intensity contribution, because of temperature-to-polarization leakage. The Commander polarization mask (CPM) is accordingly defined as the product of the (smoothed and thresholded) χ² map shown in the top panel of Fig. 41, and the low-resolution Commander CO J = 1 → 0 intensity map thresholded at 0.5 K_RJ km s^-1. The resulting mask is shown in the middle panel of Fig. 41, and excludes 17% of the sky.

The bottom panel of Fig. 41 shows a histogram of the χ² values outside the CPM83 mask, with the best-fit χ² distribution with variable number of degrees of freedom and width, fully analogous to the temperature case in Sect. 5.3. In this case, the best-fit distribution has 7.9 degrees of freedom, in excellent agreement with the theoretical expectation of 8, while the width rescaling factor that accounts for correlated noise and smoothing is 1.15, indicating that the white noise approximation underestimates the noise by 15% due to correlations.

Table 10 lists the rms of the residual maps for each frequency, analogous to Table 8 for temperature, averaged over the two Stokes Q and U parameters. The third column in this table shows the ratio between these rms values and instrumental noise; again, we observe good agreement with expectations. As for the temperature case, the values for the 30 and 353 GHz channels are significantly lower than unity, because these two frequencies dominate the synchrotron and thermal dust amplitude parameters, respectively.

Next, we assess the impact of temperature-to-polarization leakage from the CMB temperature monopole and dipole and from Galactic temperature emission by computing the synchrotron and thermal dust amplitude maps when adopting two different HFI leakage models. The first is simply the default template set adopted for the Planck 2015 release, corresponding to the results already discussed, while the second is the experimental template set described in Planck Collaboration VIII (2016) and Planck Collaboration XI (2016). From the resulting amplitude maps, we perform the following steps: first compute the polarization amplitude, P; smooth this to 3deg FWHM for synchrotron and 1deg FWHM for thermal dust; and finally compute the difference, P₂−P₁, and fractional difference, (P₂−P₁) /P₁, maps. These are shown in Fig. 42. Here we see that the absolute polarization amplitude difference between the two leakage models is around 1μK_RJ for both synchrotron and thermal dust at high Galactic latitudes, increasing to a few tens of μK_RJ in the Galactic plane. Accordingly, the fractional residuals are ≲10% at high Galactic latitudes, and they increase to about 30% in the central Galactic plane.

Finally, we comment on the polarization fractions that may be derived from these maps. First of all, we emphasize that the delivered products are maps of the Stokes Q and U parameters, not of polarization amplitude and polarization angle and fractions. This choice is primarily driven by the fact that the Stokes parameters are linear, and therefore have much simpler noise properties than the corresponding nonlinear parameters. Second, when computing polarization fractions, P/I, it is of utmost importance to recognize and account for the considerable uncertainty in this quantity with respect to the zero-level of the corresponding temperature map. To make this point concrete, we show in the top panel of Fig. 43 the naive polarization fraction derived directly from the delivered Commander thermal dust intensity and polarization maps. This map saturates the colour scale over extended regions in the southern Galactic hemisphere, nominally suggesting a polarization fraction well above 20%. However, as discussed both in Sect. 3 and Planck Collaboration VIII (2016), there is an offset in the zero-level of the zodiacal light emission of 34 μK in the current 353 GHz temperature data. Correcting for this offset results in the polarization fraction shown in the lower panel of Fig. 43, which shows significantly smaller values. Further, the raw statistical uncertainty of the 353 GHz Galactic emission zero level from H i cross-correlation alone is 23 μK (Planck Collaboration VIII 2016). The conclusion is therefore that any analysis that relies directly on the polarization fraction, as opposed to the much more stable polarization amplitude, needs to account for the significant uncertainties in the Galactic emission zero-level at 353 GHz.

6.2. Synchrotron and thermal dust angular power spectra

One of the most important goals of modern CMB cosmology is to detect primordial B-mode polarization on large angular scales, a direct observable signature of inflationary gravitational waves. The main obstacles in this search are the polarized synchrotron and thermal emission signals discussed above. In order to quantify the magnitude of this problem, we compute in this section their angular power spectra, and compare them to the expected primordial CMB spectrum. A more comprehensive analysis of the same type, but based only on the 353 GHz frequency channel, was recently published in Planck Collaboration Int. XXX (2016).

We employ the same cross-spectrum power spectrum estimator as used in Planck Collaboration Int. XXX (2016), but introduce two specific changes. First, we adopt the so-called common mask from the CMB component separation analysis presented in Planck Collaboration IX (2016), rather than the CO mask employed in the original paper, and second, we consider three different mask apodization scales (1, 2, and 5deg FWHM) as opposed to only 5deg FWHM as in Planck Collaboration Int. XXX (2016).

The EE and BB spectra resulting from the evaluation using 1deg FWHM apodization are shown in the left and right panels of Fig. 44, respectively, both plotted in terms of D_ℓ = C_ℓℓ(ℓ + 1) / 2π in thermodynamic units. Red data points show the angular power spectra for thermal dust emission at 353 GHz, and green points show synchrotron emission at 30 GHz. Each spectrum is binned with Δℓ = 25, and the plotted uncertainties are the standard deviation of the single-ℓ spectrum values within each bin. Black solid lines indicate the best-fit ΛCDM spectrum (Planck Collaboration XI 2016), and (in the BB panel only) the dashed black line shows the spectrum for a tensor-to-scalar ratio of r = 0.05. Dotted coloured lines indicate the best-fit power-law fit, D_ℓ = q (ℓ/ 80)^α, to each foreground spectrum, where the pivot scale of ℓ₀ = 80 is chosen to match that used in Planck Collaboration Int. XXX (2016). The corresponding best-fit parameters are tabulated in Table 11 for all three apodization scales, and including multipoles in the range ℓ = (10,150) for synchrotron and ℓ = (10,300) for thermal dust emission. No synchrotron results are shown for the 5deg FWHM apodization scale. In this case, the effective sky fraction is too small to allow a robust estimate of the synchrotron spectrum.

For thermal dust emission these parameters may be compared to the results presented in Planck Collaboration Int. XXX (2016), although a few caveats are in order. In particular: (1) the masks used in the two analyses are different, and the mask adopted in this paper effectively removes more sky around bright point sources after apodization; (2) the map considered in the present analysis is the Commander thermal dust map, whereas the original analysis considered the raw 353 GHz map; (3) the multipole ranges adopted for the parameter fits are slightly different; and (4) we make the fit to the single-ℓ power spectrum, not the binned spectrum.

Nevertheless, we see that the results derived here are in good agreement with those found in Planck Collaboration Int. XXX (2016). In particular, when considering the same apodization scale of 5deg FWHM, we recover an identical BB/EE ratio of 0.53 ± 0.01, and the EE power-law index agree to 0.5σ. For BB, the spectral index difference is slightly larger, but still within 2σ. The power spectrum amplitudes, on the other hand, are different because of the different effective sky fractions of the two corresponding apodized masks.

Comparing the different apodization scales, we note both that the BB/EE ratio increases slightly, and that the power-law indices steepens slightly, as the mask smoothing scale increases. This is due to thermal dust emission being a highly non-isotropic and non-Gaussian field, as discussed in Planck Collaboration Int. XXX (2016). It is not surprising that its statistical properties may vary between the Galactic plane and the high Galactic latitudes. In addition, there is an algorithmic uncertainty in the form of so-called E-to-B leakage, due to ambiguous polarization modes near the mask edge. This leakage is far stronger for foregrounds than for CMB, simply because the foreground field by construction is at its maximum near the mask boundary. As a result, it is important to specify the properties of the analysis mask when summarizing the power spectrum of a foreground field, as demonstrated in Table 11.

Overall, however, the mask dependence on the angular power spectrum is modest, and D_ℓ provides a useful summary for foreground fields as well as for the CMB field. Indeed, one of the interesting results reported by Planck Collaboration Int. XXX (2016) was the asymmetry between the B- and E-mode thermal dust power spectra, with a power ratio of BB/EE ≈ 0.5. This has strong implications for the underlying astrophysics, and indicates the presence of significant filamentary structures on intermediate angular scales. In this paper, we find that the same holds also for synchrotron emission, with an even stronger asymmetry of BB/EE ≈ 0.35. Thus, polarized synchrotron emission appears to be more strongly aligned along filamentary structures than thermal dust.

We also find similar power-law indices for synchrotron emission as for thermal dust, with α ≈ −0.4. However, the uncertainties are relatively larger, because of the lower signal-to-noise ratio of the 30 GHz channel compared to the 353 GHz channel.

These power-law fits can be used to model the total foreground level as a function of both multipole moment and frequency. This is illustrated in Fig. 45 for the 1deg FWHM apodization case in terms of iso-contour plots of the following amplitude ratio, where subscripts “s” and “d” refer to synchrotron and thermal dust. The frequency spectra, s_s(ν) and s_d(ν), are the synchrotron (GALPROP) and thermal dust (one-component greybody) spectra defined in Table 4 converted to thermodynamic units, with parameters defined by the average parameters listed in Table 5. This function is thus simply a model of the foreground-to-CMB amplitude ratio as a function of multipole and frequency.

Considering first the EE case shown in the left panel of Fig. 45, we note several interesting features. First, the horizontal ripples seen at ℓ ≳ 100 correspond to the CMB acoustic oscillations. Next, we see that the foregrounds-to-CMB ratio is smaller than unity for all multipoles above ℓ ≳ 40 for frequencies around 70 GHz, and smaller than 10% for ℓ ≳ 200. Also, recall that the corresponding power spectrum ratio goes as the square of these ratios, and we thus find that polarized foregrounds have a small effect on the EE spectrum at multipoles above a few hundred, in agreement with the results pesented in Planck Collaboration XI (2016). However, we also see that the same is by no means true at low multipoles; the foregrounds-to-CMB ratio is larger than 3 throughout the reionization peak for ℓ = 2–10.

The right panel of Fig. 45 shows the corresponding ratio for BB, but in this case two different contour sets are plotted; one for the standard ΛCDM with a vanishing tensor-to-scalar ratio (black contours), and one with a tensor-to-scalar ratio of r = 0.05 (red contours). The peak around ℓ ≈ 1000 corresponds to the signature of weak gravitational lensing, converting E-modes into B-modes, whereas the “plateau” at low multipoles in the red contours corresponds to additional primordial fluctuations from inflationary gravitational waves. First of all, we see that foregrounds are sub-dominant to the lensing signal at multipoles above ℓ ≳ 200 for frequencies around 70 GHz in this model, although they never fall below the 10% level. Second, for a vanishing tensor-to-scalar ratio the foreground-to-CMB around the recombination peak of ℓ ≈ 100 is about 3, and at the reionization peak, below ℓ ≲ 10, it is about 100. Increasing the tensor-to-scalar ratio to r = 0.05 decreases these numbers to about 2 and 20, respectively.

Before concluding this section, several caveats regarding the above observations are in order. First of all, it is important to remember that the angular power spectra reported here are computed over a large sky fraction including 72% of the sky. For a dedicated B-mode experiment, it obviously makes sense to consider more conservative masks. Second, it is also important to bear in mind that the angular spectra presented here covers only a limited multipole range, and the extrapolation to small angular scales is therefore associated with considerable uncertainty. Clearly, extrapolating actual observations that are made between ℓ ≈ 10–100 to ℓ ≈ 1000 for synchrotron emission implies strong assumptions regarding the foreground composition of both diffuse foregrounds and compact objects.

6.3. Comparison with independent data products

We now turn to consistency tests based on external (or at least independently derived) data products. Of course, given the pioneering nature of the Planck polarization observations, the number of available external cross-checks is significantly sparser compared to the temperature case. On the low-frequency side, the WMAP K-band data represents an excellent comparison for the synchrotron map, while no products of comparable data quality exists on the high-frequency side at the moment. This lack of polarized dust measurements has of course been a major limitation for the entire CMB field for a long time, and the WMAP solution to this problem was to construct a polarized dust template by combining the FDS thermal dust intensity map (Finkbeiner et al. 1999) with polarization directions from starlight polarization observations (see Page et al. 2007 for full details).

We compare our new polarized thermal dust and synchrotron maps with the WMAP maps/templates in Figs. 46 and 47, and show corresponding T–T scatter plots in Figs. 48 and 49.

Starting with the synchrotron case, we see first of all in Fig. 48 that the overall pixel-to-pixel scatter between the Planck synchrotron map and the WMAP K-band map is substantial. Indeed, based on this full-sky scatter plot, any synchrotron spectral index between β_s = −3.4 and −3.0 appears consistent with the observations. Adopting a mean value of β_s = −3.2, and assuming an effective K-band frequency of 22.6 GHz, this translates into a total scaling factor of 0.39 between K-band and 30 GHz¹⁶. This scaling factor has been applied to the K-band map shown in Fig. 46, and it also allows us to form a meaningful residual map, as seen in the bottom row of the same figure.

The relative residuals between the Planck and WMAP synchrotron maps are clearly substantial, with amplitudes reaching 5 μK at high Galactic latitudes, and with a morphology clearly associated with the scanning strategy of either Planck or WMAP, both of which have symmetries defined by the Ecliptic reference frame. Furthermore, the residuals are very large-scale in nature, and obviously dominated by the two lowest multipoles, ℓ = 2 and 3. It is therefore natural to consider what effects may cause such large-scale features. Starting with Planck, a large suite of null-tests and simulations, specifically targeting large-scale systematics in the LFI observations, is presented in Planck Collaboration II (2016). One noteworthy conclusion from that work is a significant null-test failure in the 44 GHz polarization frequency map for ℓ = 2–4, and for two 70 GHz surveys. In addition, the HFI channels between 100 and 217 GHz are also affected by low level residual systematics (Planck Collaboration VIII 2016). As a result, these observations are currently excluded from the Planck 2015 low-ℓ likelihood, which instead only relies on the 30, 70 and 353 GHz channels (Planck Collaboration XI 2016). These large-scale 44 GHz modes are likely to contribute significantly to the residuals seen in Fig. 46. For WMAP, on the other hand, the statistical uncertainties in the WMAP EEℓ = 2 and BBℓ = 3 modes are very large (Bennett et al. 2013), due to the combination of the differential detectors of WMAP, and an opening angle of 141deg between the A and B side reflectors. Although these uncertainties are appropriately described by the low-resolution WMAP covariance matrices, it is algorithmically non-trivial to account for this effect properly in component separation at higher resolution, and they are also likely to contribute to the residuals in Fig. 46.

These differences are also discussed in Planck Collaboration XXV (2016), with consistent conclusions. However, that analysis proceeds with co-adding the Planck and WMAP data sets to produce an alternative synchrotron map with higher signal-to-noise ratio, and use that map to identify new polarized synchrotron features and loops. A similar approach is not straightforward within the framework adopted in the current paper, which relies sensitively on parametric fits in the frequency domain. Within such a methodology, even just a single large-scale mode can severely bias the resulting spectral index map, rendering the full reconstruction meaningless; for a fully analogous problem, see the offset discussion in Sect. 4.1 for the temperature case. The proper way to solve this problem is either to account for the full covariance matrix in the analysis (as for instance was done by Dunkley et al. 2009), but this is only possible at very low angular resolution. An alternative approach is to marginalize over just a few modes using the template formalism described in Sect. 2. This requires accurate spatial templates in the first place, which possibly may be extracted as high-noise eigenmodes from the full noise covariance matrices. Exploratory work in this direction is already on-going, but the results have unfortunately not yet converged.

For now, we recognize that there are significant uncertainties in the synchrotron model provided, with respect to the very lowest multipoles. Fortunately, higher-ordered modes are in much better agreement between Planck and WMAP. This is important for several reasons, not least of which is estimation of the optical depth of reionization, τ, which depends sensitively on the low-ℓ modes. However, as shown in Planck Collaboration XI (2016), most of the statistical power for τ comes from ℓ = 4–6, not ℓ = 2 and 3. Thus, after removing the 44 GHz channel from the likelihood data set, all remaining null-tests pass.

A similar comparison between the Planck thermal dust map and the WMAP dust template is presented in Figs. 47 and 49. Of course, this is not a validation test of the Planck products in any way, since the two maps are by no means equivalent data products. Rather, this comparison provides an interesting reality check on the WMAP model.

The WMAP polarized dust template is provided in arbitrary units, and must therefore be re-scaled to match the amplitude of the Planck thermal dust map. The appropriate scale factor is given as the slope of the T–T correlation plot in Fig. 49. Averaging over the two Stokes parameters, we estimate this to be 480 μK_RJ per WMAP unit. This allows us to form the difference map in the bottom row of Fig. 47. With the benefit of hindsight, we see that the WMAP model is accurate to within 20–30% over most of the high-amplitude sky, although some regions show deviations at the 50% level. Larger scales are reproduced with greater fidelity than smaller scales.

Given the lack of proper external validation data sets for the polarized thermal dust emission maps, we instead compare our products with an independent internal Planck product, similar to what was done for the CO intensity maps in Sect. 5.5.3. Specifically, we compare the Commander dust map with an equivalent map derived with the SMICA algorithm (Planck Collaboration XII 2014; Planck Collaboration IX 2016). The results from this comparison are summarized in Fig. 50, in terms of: the polarization amplitude difference between the two maps (top panel); the fractional polarization amplitude difference map (middle panel); and a T–T correlation plot of the individual Stokes Q and U parameters¹⁷. First, we see that the two codes agree to better than 2% in the Galactic plane, which is very good, considering the quite different effective bandpass treatments in the two approaches. In particular, no explicit bandpass integration corrections are applied in the SMICA analysis, but all calculations are performed at the bandpass integrated frequency channel level.

However, the agreement is not equally good at high Galactic latitudes, with residuals at the roughly 1 μK_RJ level. Furthermore, these residuals have a morphology that resembles known instrumental systematics, in particular in the form of monopole and dipole leakage and ADC corrections (Planck Collaboration IX 2016). Thus, the two algorithms clearly respond differently to known systematics, and we accordingly estimate that the systematic uncertainty in the large-scale modes of the Planck polarized dust map due to instrumental effects to be (at least) 1 μK_RJ at 353 GHz. Fortunately, these instrumental effects are expected to be significantly reduced in the final Planck data release.

7. Summary and conclusions

We have presented the baseline Planck 2015 astrophysical foreground products in both temperature and polarization, as derived within a Bayesian parameter estimation framework. Combining the new Planck sky maps with complementary ancillary data in the form of the 9-yr WMAP temperature sky maps and the low-frequency 408 MHz Haslam et al. survey, we are able to reconstruct a total of six primary emission temperature mechanisms – CMB, synchrotron, free-free, spinning dust, CO, and thermal dust emission – in addition to two secondary components, namely thermal SZ emission around the Coma and Virgo clusters, and line emission between 90 and 100 GHz. For polarization, we reconstruct three primary emission mechanisms – CMB, synchrotron, and thermal dust. In addition to these astrophysical parameters, we account jointly for calibration and bandpass measurement errors, as well as monopole and dipole uncertainties. Statistical uncertainties are propagated from raw sky maps to final results by means of standard MCMC sampling techniques, while various model errors are assessed by end-to-end simulations. All data products are made publicly available, as summarized in Table 5.

Three particularly noteworthy highlights from this analysis is the following.

• We have presented the first full-sky polarized thermal dust map,which is a direct result of the exquisite sensitivity of the HFIinstrument. This map will remain a cornerstone of future CMBcosmology for the next decade or more, as the search forprimordial gravitational waves enters the next phase in whichforegrounds are more important than instrumental noise.

• We have also presented a full-sky spinning dust intensity map. In addition to its obvious scientific value, this map is also interesting for algorithmic reasons, as a clear demonstration of both the importance and power of joint global analysis. Neither WMAP nor Planck have the statistical power to disentangle spinning dust from synchrotron, but together beautiful new results emerge. We believe that this will be the default approach for virtually all future microwave surveys, as no experiment will have the power to replace Planck and WMAP by themselves. Rather, each new experiment will contribute with a new critical piece of information regarding a given phenomenon, frequency coverage, or range of angular scales, and thereby help refining the overall picture. Global Bayesian analysis provides a very natural framework for this work.

• Another useful illustration of the power of global analysis presented in this paper is the identification of important instrumental systematic errors. One example is the detection of, and correction for, systematic errors in the Planck bandpass measurements. More generally, the residual maps shown in Figs. 2, 21, and 40 comprise a treasure trove of information on instrumental systematics that should prove very valuable for improving the raw Planck sky maps before the next data release.

All things considered, the sky model presented in this paper provides an impressive fit to the current data, with temperature residuals at the few microkelvin level at high latitudes across the CMB-dominated frequencies, and with median fractional errors below 1% in the Galactic plane across the Planck frequencies. For polarization, the residuals are statistically consistent with instrumental noise at high latitudes, but limited by significant temperature-to-polarization leakage in the Galactic plane. Overall, this model represents the most accurate and complete description currently available of the astrophysical sky between 20 and 857 GHz.

Figure 51 provides an overview of the main components in both temperature (top panel) and polarization (bottom panel), summarized in terms of the brightness temperature rms evaluated over 93% and 73% of the sky, respectively. For polarization, this is the first version of such a plot that is based on observations alone. For temperature, the most recent previous version is figure 22 of Bennett et al. (2013), summarizing the WMAP temperature foreground model. While the two versions agree well in terms of total foreground power and location of the foreground minimum, there are a few subtle differences as well. The most important of these is the relative amplitude of synchrotron and spinning dust. Specifically, synchrotron dominates over spinning dust at all frequencies in the WMAP model, whereas in our new model spinning dust dominates over synchrotron between 15 and 60 GHz. Such differences are not surprising, considering the complexity of the astrophysical foregrounds at low frequencies. As emphasized repeatedly, even when combining the Planck and WMAP observations, as done in this paper, degeneracies between synchrotron, free-free and spinning dust remain the leading source of uncertainty on the low frequency side. Additional observations between, say, 2 and 20 GHz are essential to break these degeneracies. For a more complete analysis of the low-frequency foreground model presented here, we refer the interested reader to Planck Collaboration XXV (2016).

On the high-frequency side, the main outstanding issue are uncertainties in the net 545 and 857 GHz calibration, i.e., the product of calibration and bandpass uncertainties. As of today, the 545 GHz calibration is uncertain at least at the 1–2% level, and this translates into an effective 3–6% uncertainty for the 857 GHz channel in our fits (in order to maintain a physical thermal dust frequency scaling). Cross-correlations with H i observations suggests a total systematic error on the thermal dust temperature at high Galactic latitudes of 1–2 K. Recognizing both calibration and modelling errors, we emphasize that the thermal dust model presented here does not represent an accurate model of frequencies beyond 857 GHz. For instance, naive extrapolation to the DIRBE 100 μm channel results in residuals as large as 40%. Both more physical models and better calibration are needed to extend into this regime. In addition, it is important to note that the current model makes no attempt at separating Galactic thermal dust emission from CIB fluctuations, and these therefore constitute a significant contaminant in our thermal dust model on small angular scales.

Finally, for polarization the main limitations are instrumental systematics, primarily in the form of temperature-to-polarization leakage, uncertainties in the analogue-to-digital conversion, and very long time constants (Planck Collaboration I 2016; Planck Collaboration VI 2016; Planck Collaboration VIII 2016). Thus, although the new Planck 2015 observations already have opened up a completely new window on the physics of our own Galaxy, through its deep observations of polarized thermal dust, more work is required in order to fully realize the science potential of the Planck measurements. This will be the main focus of the Planck analysis efforts in the coming months.

---

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration

References

All Tables

All Figures