© ESO, 2014

1. Introduction

This paper, one of a set associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration I 2014), presents a new parametrization of dust emission that covers the whole sky, at 5′ resolution, based on data from 353 to 3000 GHz (100 to 850 μm).

Because it is well mixed with the gas and because of its direct reaction to UV photons from stars, dust is a great tracer of the interstellar medium (ISM) and of star formation activity. On the other hand, for many studies in extragalactic astrophysics and cosmology, Galactic interstellar dust is a nuisance, a source of extinction and reddening for UV to near-infrared observations and a contaminating emission in the infrared to millimetre wavelengths. Thanks to the sensitivity, spectral coverage, and angular resolution of Planck, this model of dust emission brings new constraints on the dust spectral energy distribution (SED), on its variations across the sky, and on the relationships between dust emission, dust extinction, and gas column density. In particular, this model of dust emission provides a new map of dust extinction at 5′ resolution, aimed at helping extragalactic studies.

The emission in the submillimetre range arises from the bigger dust grains that are in thermal equilibrium with the ambient radiation field. Thermal dust emission is influenced by a combination of the dust column density, radiation field strength, and dust properties (size distribution, chemical composition, and the grain structure). When the effect of the radiation field can be estimated (using the dust temperature as a probe) and the dust properties assumed, the dust optical depth is possibly the most reliable tracer of interstellar column density, and therefore of mass for objects at known distances. Dust optical depth is used to estimate the mass of interstellar clumps and cores (Ossenkopf & Henning 1994) in particular with the higher resolution Herschel data (Launhardt et al. 2013), to study the statistical properties of the ISM structure and its link with gravity, interstellar turbulence, and stellar feedback (Peretto et al. 2012; Kainulainen et al. 2013), and as a way to sample the mass of the ISM in general (Planck Collaboration XIX 2011). The accuracy of these determinations depends on, among other things, the frequency range over which the dust spectrum is observed. The combination of Planck and IRAS data offers a new view on interstellar dust by allowing us to sample the dust spectrum from the Wien to the Rayleigh-Jeans sides, at 5′ resolution over the whole sky.

Dust emission, with extinction and polarization, is a key element to constrain the properties of interstellar dust (Draine & Li 2007; Compiègne et al. 2011). The dust emissivity (i.e., the amount of emission per unit of gas column density) and the shape of the dust SED provide information on the nature of the dust particles, in particular their structure, composition, and abundance, related to the dust-to-gas ratio.

Changes in the dust emissivity and the shape of the dust SED can be related to dust evolutionary processes. Interstellar dust grains are thought to be the seeds from which larger particles form in the ISM, up to planetesimals in circumstellar environments (Brauer et al. 2008; Beckwith et al. 2000; Birnstiel et al. 2012). This growth of solids can be followed in earlier phases of the star-formation process, at the protostellar phase and even before, in molecular clouds and in the diffuse ISM. Many studies have revealed increases of the dust emissivity with (column) density in molecular clouds accompanied by a decrease in dust temperature (Stepnik et al. 2003; Schnee et al. 2008; Planck Collaboration XXV 2011; Arab et al. 2012; Martin et al. 2012; Roy et al. 2013). One explanation is that grain structure is changing through aggregation of smaller particles, enhancing the opacity (Köhler et al. 2011). Planck’s spectral coverage allows us to model the big grain thermal emission, in particular its spectral index that is related to the grain composition and structure (Ormel et al. 2011; Meny et al. 2007; Köhler et al. 2012). Because of its full-sky coverage, Planck can also reveal variations of the dust SED with environment, enabling us to better understand the evolutionary track of dust grains through the ISM phases.

Dust emission is one of the major foregrounds hampering the study of the cosmic microwave background (CMB). The thermal dust emission peaks at a frequency close to 2000 GHz but its emission is still a fair fraction of the CMB anisotropies in the 20–200 GHz range where they are measured. This is even more the case in polarization (Miville-Deschênes 2011). The model of dust emission proposed by Finkbeiner et al. (1999) based on data from previous satellite missions (IRAS and COBE) made an important contribution to the field, in guiding the design of CMB experiments and in helping the data analysis by providing a spatial template and a spectral dependence of the dust emission at CMB frequencies. It is still the basis of recent models of Galactic foreground emission (Delabrouille et al. 2013). With its frequency coverage that bridges the gap between IRAS and the CMB range, its high sensitivity, and its better angular resolution, Planck offers the opportunity to develop a new model of thermal dust emission.

Estimating reddening and extinction by foreground interstellar dust is a major issue for observations of extragalactic objects in the UV to near-infrared range. Major efforts have been made toward producing sky maps that provide a way to correct for the chromatic extinction of light by Galactic interstellar dust on any line of sight. First Burstein & Heiles (1978) used H i as a proxy for dust extinction by correlating integrated 21 cm line emission with extinction estimated from galaxy counts. It was subsequently discovered that H i is not a reliable tracer of total column density N_H for N_{H i} greater than a few 10²⁰ cm^-2 due to molecular gas contributions (Lebrun et al. 1982; Boulanger & Pérault 1988; Désert et al. 1988; Heiles et al. 1988; Blitz et al. 1990; Reach et al. 1994; Boulanger et al. 1996). It was then proposed to use dust emission as a more direct way to estimate dust extinction. By combining 100, 140, and 240 μm data (DIRBE and IRAS) Schlegel et al. (1998) produced an all-sky map of dust optical depth at 100 μm that was then calibrated into dust reddening by correlating with colour excesses measured for galaxies. The work presented here is the direct continuation of these studies. Like Schlegel et al. (1998) we also propose a map of E(B − V) based on a model of dust emission calibrated using colour excess measurements of extragalactic objects, here quasars.

The paper is organized as follows. The data used and the preprocessing steps are presented in Sect. 2. The model of the dust emission, SED fit methodology, the exploration of potential biases, and the all-sky maps of dust parameters are described in Sect. 3. The results of the Galactic dust model are analysed in Sect. 4. Sections 5 and 6 describe specifically how the dust emission model compares to other tracers of column density. The Planck dust products, the dust model maps, and the E(B − V) map aimed at helping extragalactic studies to estimate Galactic extinction are detailed in Sect. 7 and compared with similar products in the literature. Concluding remarks are given in Sect. 8.

2. Data and preprocessing

The analysis presented here relies on the combination of the Planck data from the HFI instrument at 857, 545, and 353 GHz (respectively 350, 550, and 850 μm) with the IRAS 100 μm (3000 GHz) data.

2.1. Planck data

For Planck we used the HFI 2013 delivery maps (Planck Collaboration VI 2014), corrected for zodiacal emission (ZE – see Planck Collaboration XIV 2014). Each map was smoothed to a common resolution of 5′, assuming a Gaussian beam². The 353 GHz map, natively built in units of K_CMB, was transformed to MJy sr^-1 using the conversion factor given by Planck Collaboration IX (2014). The CMB anisotropies map provided by the SMICA algorithm (Planck Collaboration XII 2014), which has an angular resolution of 5′, was removed from each Planck HFI map.

As shown in Planck Collaboration XIII (2014), ¹²CO and ¹³CO rotational lines fall in each of the Planck HFI filters, except at 143 GHz. At 857 and 545 GHz the CO lines (J = 5 → 4 and J = 4 → 3, respectively) are very faint compared to the dust emission and they are not considered here. On the other hand, emission from the ¹²CO J = 3 → 2 line was detected in the 353 GHz band (Planck Collaboration XIII 2014). Nevertheless, this emission is still faint compared to the dust emission whereas the noise on the Planck CO emission estimate in the 353 GHz band is quite high (see Planck Collaboration XIII 2014 for details). The detection of the ¹²CO J = 3 → 2 line emission by Planck is above 3σ for only 2.6% of the sky. When detected above 3σ, this emission is on average 2% of the 353 GHz specific intensity. It contributes 5% or more of the 353 GHz specific intensity for only 0.3% of the sky. Given such a relatively small contribution we did not subtract CO emission from the data so as not to compromise the 353 GHz map through the adverse impact of the noise of the ¹²CO J = 3 → 2 product.

2.2. The 100 μm map

The 100 μm map used in this analysis is a combination of the IRIS map (Miville-Deschênes & Lagache 2005) and the map of Schlegel et al. (1998, hereafter SFD), both projected on the HEALPix³ grid (Górski et al. 2005) at N_side = 2048. Both IRIS and SFD maps were built by combining IRAS and DIRBE 100 μm data. Nevertheless these two maps show differences at large scales due to the different assumptions used for the ZE removal. Miville-Deschênes & Lagache (2005) used the DIRBE 100 μm map from which ZE was removed by the DIRBE team, using the model of Kelsall et al. (1998). On the other hand, Schlegel et al. (1998) used their own empirical approach to remove ZE based on a scaling of the DIRBE 25 μm data. Because it is based on data and not on a model, the SFD correction is closer to the complex structure of the ZE and provides a better result. This can be assessed by looking at the correlation of the IRIS and SFD maps with H i in the diffuse areas of the sky (1 <N_{H i}< 2 × 10²⁰ cm^-2), as detailed in Appendix A.1. The uncertainty of the slope of the correlation with N_{H i} and the standard deviation of the residual is about 30% lower for the SFD map compared to the IRIS map. For that reason (and others described in Appendix A.1) we favour the use of the SFD map at large scales.

At scales smaller than 30′, the IRIS map has several advantages over the SFD map⁴. The IRIS map is at the original angular resolution (4.́3) of the IRAS data while SFD smoothed the map to 6.́1. IRIS also benefits from a non-linear gain correction that is coherent for point sources and diffuse emission. Finally point sources were kept in the IRIS map whereas SFD removed some of them (mostly galaxies but also ISM clumps). To combine the advantages of the two maps, we built a 100 μm map, I₁₀₀, that is compatible with SFD at scales larger than 30′ and compatible with IRIS at smaller scales: ${\mathit{I}}_{\mathrm{100}}\mathrm{=}{\mathit{I}}_{\mathrm{IRIS}}\mathrm{-}{\mathit{I}}_{\mathrm{IRIS}}\mathrm{\otimes }{\mathit{f}}_{\mathrm{IRIS}}^{\mathrm{30}}\mathrm{+}{\mathit{I}}_{\mathrm{SFD}}\mathrm{\otimes }{\mathit{f}}_{\mathrm{SFD}}^{\mathrm{30}}\hspace{0.17em}\mathit{,}$(1)where I_IRIS and I_SFD are, respectively, the IRIS and the SFD maps, and ${\mathit{f}}_{\mathit{i}}^{\mathrm{30}}$ is the complementary Gaussian kernel needed to bring the maps to 30′ resolution.

2.3. Zero level

The fit of the dust emission requires that the specific intensity at each frequency and at each sky position is free of any other emission. In particular the zero level of each map should be set in such a way that it contains only Galactic dust emission. In order to set the zero level of the maps to a meaningful Galactic reference we applied a method based on a correlation with H i, as described in Planck Collaboration VIII (2014).

Some precautions need to be taken here because the ratio of dust to H i emission might vary locally due to variations of the radiation field or of the dust optical properties. Locally the 21 cm emission might not be a perfect tracer of the column density due to H i self-absorption effects or to the presence of ionized or molecular gas. Nevertheless, the correlation between dust and H i emission is known to be tight in the diffuse ISM where most of the gas is atomic. This correlation has been used several times to establish the dust SED (Boulanger et al. 1996; Planck Collaboration XXIV 2011), to isolate the cosmic infrared background (CIB; Puget et al. 1996; Planck Collaboration XVIII 2011; Pénin et al. 2012), and to establish a Galactic reference for dust maps (Burstein & Heiles 1978; Schlegel et al. 1998).

The excess of dust emission with respect to the H i correlation has been used to reveal gas in molecular form, even in regions where CO emission was not detected (Désert et al. 1988; Blitz et al. 1990; Reach et al. 1998; Planck Collaboration XIX 2011). Such an excess can be observed at column densities as low as N_{H i} = 2 × 10²⁰ cm^-2. Using this as an upper limit for our correlation studies also ensures that self-absorption in the 21 cm line emission is not important. Note that this is also below the threshold at which significant H₂ is seen in the diffuse ISM (Gillmon et al. 2006; Wakker 2006; Rachford et al. 2002, 2009).

To estimate the Galactic reference of the IRAS and Planck data, the maps were correlated against the 21 cm LAB data (Kalberla et al. 2005), integrated in velocity. The ranges are referred to as LVC, low velocity gas with | v_LSR | < 35 km s^-1, and IVC, intermediate velocity gas with 35 < | v_LSR | < 70 km s^-1 (Albert & Danly 2004). HVC, high-velocity clouds with | v_LSR | > 70 km s^-1 are excluded. For LVC and IVC separately, column density maps N_{H i} assuming optically thin emission are given in Fig. 1 in the all-sky Mollweide equal-area projection and in Fig. 2 in a complementary polar orthographic projection⁵.

For the correlation, all data sets were convolved to 1° resolution and projected on an N_side = 128 grid. The correlation was performed using only pixels with the LVC N_{H i}< 2 × 10²⁰ cm^-2, discarding pixels with detected IVC above 0.1 × 10²⁰ cm^-2. The resulting area covers 11.5% of the sky, corresponding to more than 4700 deg². This mask is presented in both the all-sky view in Fig. 3 (left panel) and the polar view in Fig. 4. We refer to this throughout as the “low N_{H i} mask”.

The correlations are shown in Fig. 5, left. The values of the offset that were removed from the 3000 and 857 GHz maps are given in Table 1. We have checked that the offsets estimated in that way are not sensitive to the resolution of the H i data or to the area of the sky selected. We have also checked that the assumption that the 21 cm emission is optically thin does not introduce a significant bias in the analysis. For example, on assuming T_spin = 80 K in converting 21 cm emission to N_{H i} (Lockman & Condon 2005), the changes in the offsets at 857 and 3000 GHz are only about 0.02 MJy sr^-1. This is as expected because in the LAB data for the diffuse areas of the sky considered here values of the 21 cm line brightness temperature higher than 10 K are exceptional. Compatible offsets, within the quoted uncertainties, were found using 16′ Galactic All Sky Survey 21 cm data (McClure-Griffiths et al. 2009) of the area around the Galactic south pole (Planck Collaboration Int. XVII 2014) and 9′ data obtained on smaller regions in the northern sky at the Green Bank Telescope (Planck Collaboration XXIV 2011). Finally, in Fig. 5 we note a systematic excess of the dust emission at 857 and 3000 GHz with respect to the correlation at the lowest N_{H i}. This is also seen at 545 and 353 GHz and it is discussed further in Sect. 5.4.

The correlation of dust emission between Planck frequencies is observed to be tight (the correlation of Planck is less tight with the IRAS 3000 GHz map). We took advantage of this to estimate more precisely the Galactic zero level of the 353 and 545 GHz channels. They were obtained by correlation with the 857 GHz map on a larger mask with total LVC plus IVC N_{H i}< 3 × 10²⁰ cm^-2 (Fig. 3, right panel). These correlations are shown in Fig. 5 (right). The offset values obtained in this way (see Table 1) are compatible within 1σ with the offset values deduced from the H i correlation in the smaller mask (0.104 and 0.088 MJy sr^-1 respectively at 545 and 353 GHz). At these frequencies we favour the offset obtained with the correlation with 857 GHz as it minimizes the aforementioned effects in the correlation between dust and gas emission.

In the process, faint dipole residuals were identified in the 353 and 545 GHz maps. The orientation of these residual dipoles coincides with the solar dipole and their amplitudes (see Table 1) corresponds to + 7.6% and − 0.9% of the solar dipole amplitude at 545 and 353 GHz respectively, which is within the calibration uncertainties at these frequencies. They were removed prior to the dust SED fit.

3. Model of the dust emission

3.1. Dust emission observed by Planck

The emission from interstellar dust in the far-infrared (FIR) to millimetre range is dominated by the emission from the biggest grains that are in thermal equilibrium with the local radiation field. Many studies and reviews have been dedicated to this subject (e.g., Draine 2003; Draine & Li 2007; Compiègne et al. 2011).

In the optically thin limit, the SED of emission from a uniform population of grains is well described, empirically, by a modified blackbody (MBB): ${\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathit{\tau }}_{\mathit{\nu }}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\hspace{0.17em}\mathit{,}$(2)where I_ν is the specific intensity, B_ν(T) is the Planck function for dust at temperature T, and τ_ν is the frequency-dependent dust optical depth modifying the blackbody shape of the SED. The optical depth is the product of the dust opacity, σ_{e ν}, and the gas column density, N_H: ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{\sigma }}_{\mathrm{e}\hspace{0.17em}\mathit{\nu }}\hspace{0.17em}{\mathit{N}}_{\mathrm{H}}\hspace{0.17em}\mathit{.}$(3)Alternatively, the optical depth is the product of the dust emissivity cross section per unit mass, κ_ν (in cm² g^-1), and the dust mass column density, M_dust: ${\mathit{\tau }}_{\mathit{\nu }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}\hspace{0.17em}{\mathit{M}}_{\mathrm{dust}}\hspace{0.17em}\mathit{,}$(4)where M_dust = rμm_HN_H, with r being the dust-to-gas mass ratio, μ the mean molecular weight, and m_H the mass of a hydrogen atom. Note that κ_ν depends on the chemical composition and structure of dust grains, but not the size for particles small compared to the wavelength, as here. It is usually described as a power law κ_ν = κ₀(ν/ν₀)^β (Hildebrand 1983; Compiègne et al. 2011), where κ₀ is the emission cross-section at a reference frequency ν₀. Put together, the emission of dust of a given composition and structure and in thermal equilibrium is: ${\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathit{\kappa }}_{\mathrm{0}}\hspace{0.17em}\mathrm{(}\mathit{\nu }\mathit{/}{\mathit{\nu }}_{\mathrm{0}}{\mathrm{)}}^{\mathit{\beta }}\hspace{0.17em}\mathit{r}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}{\mathit{m}}_{\mathrm{H}}\hspace{0.17em}{\mathit{N}}_{\mathrm{H}}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{(}\mathit{T}\mathrm{)}\hspace{0.17em}\mathit{.}$(5)In practice the shape of the observed SED depends on three main parameters. First, the equilibrium temperature is set by the radiation field strength, parametrized by the scaling factor, U, of the mean interstellar radiation field (ISRF) in the solar neighbourhood from Mathis et al. (1983); note that in dense regions, U is decreased because of attenuation. Second, the grain size distribution (Mathis et al. 1977; Weingartner & Draine 2001) is important; exposed to the same ISRF, bigger grains have a lower equilibrium temperature than smaller ones. Third, the dust structure and composition determine not only the optical and UV absorption cross section, but also the emission cross-section, the frequency-dependent efficiency to emit radiation, usually modelled as above as a power law (κ₀ν^β) but possibly more complex depending on dust properties (β could vary with frequency and/or grain size and/or grain temperature). In a given volume element along the line of sight, the distribution of dust grain sizes will naturally create a distribution of equilibrium temperatures. In addition, dust properties might vary along the line of sight. Furthermore, U might also change along some lines of sight. Therefore, the observed dust SED is a mixture of emission modified by these effects, the sum of several different MBBs. Nevertheless, the simplification of fitting a single MBB is often adopted and indeed here, with only four photometric bands available, is unavoidable. The parametrization of the MBB for the empirical fit is: ${\mathit{I}}_{\mathit{\nu }}\mathrm{=}{\mathit{\tau }}_{{\mathit{\nu }}_{\mathrm{0}}}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{obs}}\mathrm{\right)}\hspace{0.17em}{\left(\frac{\mathit{\nu }}{{\mathit{\nu }}_{\mathrm{0}}}\right)}^{{\mathit{\beta }}_{\mathrm{obs}}}\hspace{0.17em}\mathit{,}$(6)where ν₀ is a reference frequency at which the optical depth τ_ν0 is estimated (ν₀ = 353 GHz in our SED applications in this paper).

The main challenge is then to relate the parameters of the fit to physical quantities. It has been shown by many authors (Blain et al. 2003; Schnee et al. 2007; Shetty et al. 2009b; Kelly et al. 2012; Juvela & Ysard 2012a,b; Ysard et al. 2012) that, in general, the values of T_obs and β_obs recovered from an MBB fit cannot be related simply to the mass-weighted average along the line of sight of the dust temperature and spectral index. Even for dust with a spectral index constant in frequency (i.e., β does not depend on ν), the distribution of grain sizes and the variations of U along the line of sight could introduce a broadening of the SED relative to the case of a single dust size and single U. In addition, the dust luminosity is proportional to T^{4 + β} and so dust that is hotter for any reason, including efficiency of absorption, will contribute more to the emission at all frequencies than colder dust. Therefore, the observed SED is not a quantity weighted by mass alone. The dust SED is wider than a single MBB due to the distribution of T, and so the fit is bound to find a solution where β_obs<β and, in consequence, where T_obs is biased toward higher values. This results in dust optical depth that is generally underestimated: τ_obs<τ. This effect is somewhat mitigated when lower frequency data are included. In the Rayleigh-Jeans limit the effect of temperature is low and the shape of the spectrum is dominated by the true β. For T = 15–25 K dust, this range is at frequencies lower than ν = kT/h = 310–520 GHz.

Models like the ones of Draine & Li (2007) and Compiègne et al. (2011) go beyond the simple MBB parametrization by incorporating the variation of the equilibrium temperature of grains due to the size distribution. The model of Draine & Li (2007) also includes a prescription for the variation of U along the line of sight, but assumes fixed dust properties. Nevertheless, there are still many uncertainties in the properties of dust (the exact size distribution of big grains, the optical properties, and the structure of grains), in the evolution of these properties from diffuse to denser clouds, and in the variation of the radiation field strength along the line of sight.

Therefore, for our early exploration of the dust SED over the whole sky, at 5′ resolution and down to 353 GHz (850 μm), we believe that it is useful to fit the dust SED using the empirical MBB approach, before attempting to use more physical models that rest on specific hypotheses. The three parameters τ_ν0, T_obs and β_obs obtained from the MBB fit should be regarded as a way to fit the data empirically; the complex relationship between these recovered parameters and physical quantities needs to be investigated in detail with dedicated simulations (e.g., Ysard et al. 2012), but is beyond the scope of this paper.

3.2. Implementation of the SED fit

The fit of the dust SED with a MBB model has been carried out traditionally using a χ² minimization approach. Recently, alternative methods for fitting observational data with limited spectral coverage have been proposed, based on Bayesian or hierarchical models (Veneziani et al. 2010, 2013; Kelly et al. 2012; Juvela et al. 2013). These new methods were developed specifically to limit the impact of instrumental noise on the estimated parameters. Even though these methods offer interesting avenues, we developed our own strategy to fit the dust SED over the whole sky because of another challenge to be mitigated, arising from the cosmic infrared background anisotropies (the CIBA). Although this has been overlooked, it can be dominant in the faint diffuse areas of the sky, as we demonstrate. We proceeded with a method based on the standard χ² minimization (see Appendix B) but implemented a two-step approach that limits the fluctuations of the estimated parameters at small angular scales induced by noise and the CIBA. In developing the methodology we have explored using data degraded to lower resolution and smaller N_side.

3.2.1. Frequency coverage

One possible source of bias in the fit is the number of bands and their central frequency. The combination of Planck 353 to 857 GHz and IRIS 3000 GHz data allows us to sample the low and high frequency sides of the dust SED. For a typical temperature of 20 K, the peak of the emission is at a frequency of 2070 GHz. This falls in a gap in the frequency coverage, between 857 GHz and 3000 GHz. It is thus a concern that a fit of the Planck and IRIS data might bias the recovered parameters T_obs and β_obs. To explore this we combined the Planck data with the DIRBE data at 1250, 2143, and 3000 GHz (100, 140, and 240 μm), all smoothed to 60′, providing a better sample of the dust SED near its peak. We found that the recovered dust parameters are stable whether DIRBE data are used or not; no bias is observed in T_obs, β_obs, and τ₃₅₃ compared with results obtained using just Planck and IRIS data.

We also evaluated the potential advantage of fitting the SED with only the Planck 353 to 857 GHz data, a more coherent dataset not relying on the IRIS data. However, because the Wien part of the SED is not sampled the results showed a clear bias of T_obs, toward lower values; consequently, when extrapolated to 100 μm, the fits greatly underestimate the emission detected in the IRIS data. Therefore, in the following the χ² minimization fit was carried out on the data described in Sect. 2: the 857, 545, and 353 GHz Planck maps, corrected for zodiacal emission, and the new 100 μm map obtained by combining the IRIS and SFD maps.

3.2.2. Noise and cosmic infrared background anisotropies

Degeneracy (anticorrelation) of the estimated T_obs and β_obs, inherent to the MBB fit of dust emission in the presence of noise, has had dedicated specific study (Shetty et al. 2009a; Juvela & Ysard 2012a). As mentioned above, the CIBA is also a contaminating source in the estimate of the MBB parameters.

The CIBA is the result of the combined emission of distant unresolved galaxies. Its structure on the sky is diffuse. The angular power spectrum, with C_ℓ ∝ ℓ^-1 for 100 <ℓ< 2000 according to Planck Collaboration XVIII (2011) and Planck Collaboration XXX (2014), reveals the large-scale structure of the Universe at high redshift. The zero levels of the maps were set through correlation with H i and so the data used in our study are insensitive to the monopole of the CIBA. However, the anisotropies, or fluctuations, are present in the maps. Because the CIBA power spectrum is flatter than that of interstellar dust emission (C_ℓ ∝ ℓ^-2.9, Miville-Deschênes et al. 2007), in relative terms the CIBA is more visible at small scales. Another feature of the CIBA is that its structure on the sky is correlated in frequency, though only partially because galaxies at different redshifts contribute to the emission at different frequencies. Because of this partial correlation in frequency, the CIBA cannot be treated in the same way as instrumental noise in the fit. But it cannot be included as another component in the fitting function either. Nevertheless, the CIBA has an impact on the parameters of the fit; like the instrumental noise, the CIBA introduces an anticorrelation between T_obs and β_obs.

One option to limit the effect of noise and the CIBA is to reduce the number of free parameters in the fit. In that context we have examined the possibility of fitting the dust SED over the whole sky, at 5′ resolution, using a fixed β_obs with values between 1.5 and 1.8. A value of β_obs = 1.65 provides the best fit with a reduced χ² lower than unity everywhere on the sky but this is mostly due to the fact that we took into account calibration uncertainties in the fit (Appendix B). What is statistically significant is the fact that on about 25% of the sky the reduced χ² is improved by letting β_obs be a free parameter. This happens mostly in bright regions of the sky where the noise is not an issue. In molecular clouds and in the Galactic plane, there are variations in the shape of the SED that cannot be fit with only two parameters.

While fixing a parameter of the fit over the whole sky might be too strict, it might not be necessary to have all three parameters at full resolution to describe the data. We have thus evaluated the possibility of estimating one of the parameters at a lower resolution than the others. In the following we explore the impact of noise and the CIBA on the parameters of the fit as a function of angular resolution.

Both the noise and the CIBA have flatter power spectra than dust emission, and so we expect the intrinsic dust parameters T_obs, β_obs, and τ₃₅₃ to have a smoother structure on the sky than noise and the CIBA, except perhaps in bright photon-dominated regions where the shape of the dust SED might vary rapidly at small scales due to radiative transfer effects and potentially fast dust evolution. Smoothing the maps by different amounts before fitting on each pixel therefore offers the advantage of revealing both this spatially smoother solution and the important impact of noise and the CIBA on the result of the fit. This is illustrated in Fig. 6 where we present normalized distribution functions (NDFs) of T_obs, β_obs, and τ₃₅₃ obtained with data smoothed to 5′, 15′, 30′, 60′, and 120′, selecting only pixels corresponding to the low N_{H i} mask (Fig. 3) to highlight a regime of relatively low signal-to-noise ratio. Smoothing the data has no real impact on the average value of the parameters, but the standard deviations of T_obs and β_obs go down rapidly with smoothing; at 5′ resolution the standard deviation of T_obs and β_obs is about twice as large as with smoothed data. On the other hand, the dispersion of τ₃₅₃ is less affected by smoothing; it is dominated instead by cosmic variance, the considerable range of column densities even within this low N_{H i} mask⁶. Of course the dust parameters might also vary at small scales and so a trade-off needs to be found.

To explore and quantify the impact of noise and the CIBA on the fit at different angular resolutions, for later comparison with the actual dispersion, we used Monte Carlo simulations of the SED, including dust emission, noise, and the CIBA. The details of the Monte Carlo simulations, including the information on the inter-frequency coherence, are described in Appendix C. We considered five different angular resolutions of the data: 5′, 15′, 30′, 60′, and 120′. The noise and CIBA levels used for each resolution are given in Table C.1. Here we present results for an SED appropriate to Fig. 6 by adopting the median dust parameters found in the low N_{H i} mask that corresponds to the faintest 10% of the sky. We simulated 10⁵ realizations of this SED to which noise and the CIBA were added. For each realization the three parameters τ₃₅₃, T_obs, and β_obs were estimated as in Appendix B. The 1σ dispersions of T_obs, and β_obs obtained at each resolution are given in Table 2, for noise and the CIBA separately. The simulated effect of smoothing on the β_obs – T_obs anticorrelation is shown in Fig. 7.

At full resolution the noise is the dominant source of error on the retrieved parameters. For the specific faint dust spectrum considered here, the noise produces an uncertainty of 2.1 K while the uncertainty due to the CIBA is only 0.39 K. The same is true for β_obs: the uncertainties are 0.49 and 0.11 for the noise and the CIBA, respectively. However, even with moderate smoothing of the data, the impact of noise on the fit reduces sharply, whereas the reduction of the impact of the CIBA is less dramatic. This arises because the CIBA has a power spectrum that is steeper than that of typical (white) noise. In addition, unlike the CIBA, noise has power up to the pixel scale (i.e., it is not attenuated by the beam). For example, as seen in Table 2, for data smoothed to 30′, the noise levels of the Planck and IRAS data go down by a factor 18.5 while the CIBA standard deviation decreases only by a factor 2.5. As a result, our simulations show that for data smoothed to resolution larger than 15′ the CIBA becomes the main source of error.

Similar relative effects as a function of resolution are seen in the β_obs – T_obs anticorrelation in Fig. 7. At all resolutions the estimates of T_obs and β_obs lie within an ellipse in β_obs – T_obs space. The orientation and extent of the ellipse depends on the amplitudes of the noise and of the CIBA, which are both different and in a different ratio at each resolution. In all cases the ellipse is centred on the input values, demonstrating that noise and the CIBA do not bias the estimate of T_obs and β_obs. This is the case even though the CIBA has a flatter (broader) SED than interstellar dust. Because the CIB monopole was removed in the data and therefore not included in the simulation, the CIBA produces as many negative as positive CIB fluctuations on the sky at each frequency. Because they are (partially) correlated in frequency, positive CIB fluctuations bias the SED and descriptive dust parameters toward a flatter SED while negative CIB fluctuations have the opposite effect, toward a steeper SED.

3.2.3. The two-step approach

Given the impact of noise and the CIBA on the recovered parameters, described in the previous section, we have chosen to fit the data in two steps. First, we fit the data smoothed to 30′ (but on the N_side = 2048 grid). As shown in Figs. 6 and 7 this greatly reduces the effect of noise on the estimate of β_obs and T_obs. Second, we fit the data at 5′ resolution with a fixed β_obs taken from the map of β_obs obtained with data at 30′ resolution. That way two degrees of freedom (τ₃₅₃ and T_obs – see Eq. (6)) are still available to capture the variations of the dust SED at full resolution while limiting the effect of the β − T degeneracy due to noise.

This two-step approach is in the same spirit as the one implemented in the Commander-Ruler algorithm (Planck Collaboration XII 2014). The advantage of such methods arises by favouring a spatially smoother solution for parameters that are not expected to vary strongly at small scale. In the second fit we chose to fix β_obs rather than T_obs. It is not yet clear how the actual spectral index of the grain opacity, β, might vary on small scales (some models even assume that it is constant: Draine & Li 2007; Compiègne et al. 2011). However, the dust temperature is expected to vary on small scales, especially in dense regions of the ISM due to the attenuation of the radiation field.

We performed Monte Carlo simulations to evaluate the contributions of noise and the CIBA (partly correlated in frequency – see Appendix C) to variations, whence uncertainties, of the recovered T_obs and β_obs for the specific case of the adopted two-step fit (30′ and 5′). Figure 8 illustrates the uncertainties of T_obs and β_obs arising from noise and the CIBA for dust SEDs on lines of sight with increasing τ₃₅₃.

For the typical SED (τ₃₅₃ = 9.6 × 10^-7, T_obs = 20.8 K, and β_obs = 1.55) corresponding to the 10% faintest area of the sky, we made a comparison of the Monte-Carlo results for the two-step fit and the direct three-parameter fit (Sect. 3.2.2). The uncertainties of the direct fit at 5′ are δ(T_obs) = 2.1 K and δ(β_obs) = 0.50, adding the contributions of noise and CIBA in quadrature (see Table 2). For the same SED parameters, the uncertainties of the two-step fit are δ(T_obs) = 0.8 K and δ(β_obs) = 0.06 (see Fig. 8). In addition, the results of the two-step fit simulations indicate that for both T_obs at 5′ and β_obs at 30′ the CIBA has a greater contribution than the noise, contrary to the situation for the direct three-parameter fit.

For the faintest 0.4% of the sky (top axis in Fig. 8), the results of the simulations indicate that the combined effects of noise and the CIBA produce variations δT_obs = 3.0 K and δβ_obs = 0.2. On the other hand, for about 93% of the sky the variations are much smaller, δT_obs< 1.0 K and δβ_obs< 0.1, i.e., <5% and <6% fractional error, respectively. This is in accordance with the fact that at 353 GHz, where the CIBA is the strongest contaminant, about 93% of the sky has I₃₅₃> 3σ_CIBA(353).

3.3. Parameters and uncertainties

The all-sky maps of the dust parameters, T_obs, β_obs, and τ₃₅₃, and of their fractional uncertainties are presented in Figs. 9 and 10, respectively. The precision of the three parameters is of the order of a few percent on most of the sky. The uncertainties shown here are based on the statistical ones returned by the χ² minimization fit assuming that the model is a good representation of the data. For β_obs the uncertainty is from the 30′ fit. For T_obs we added quadratically the fractional uncertainties from the 30′ and 5′ fits to include the covariance between T_obs and β_obs, whence $\mathit{\delta }{\mathit{T}}_{\mathrm{obs}}\mathrm{=}{\mathit{T}}_{\mathrm{obs}}\hspace{0.17em}\sqrt{{\left(\frac{\mathit{\delta }{\mathit{T}}_{\mathrm{obs}\mathit{,}\mathrm{30}}}{{\mathit{T}}_{\mathrm{obs}\mathit{,}\mathrm{30}}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{\delta }{\mathit{T}}_{\mathrm{obs}\mathit{,}\mathrm{5}}}{{\mathit{T}}_{\mathrm{obs}\mathit{,}\mathrm{5}}}\right)}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}$(7)where the subscripts 5 and 30 refer to the parameter or uncertainty maps obtained at 5′ and 30′, respectively. Similarly, the uncertainty of τ₃₅₃ is from the quadratic sum of the fractional uncertainties of ${\mathit{I}}_{\mathrm{353}}^{\mathrm{m}}$ and B₃₅₃(T_obs) where ${\mathit{I}}_{\mathrm{353}}^{\mathrm{m}}$ is the reconstructed model of the emission at 353 GHz. To estimate the uncertainty of B₃₅₃(T_obs), we simply computed max | B₃₅₃(T_obs ± δT_obs) − B₃₅₃(T_obs) |.

The three uncertainty maps have a similar spatial structure. In general the fractional uncertainties are higher in the most diffuse areas of the sky (where the noise and the CIBA have a more important contribution) and in the inner Galaxy region. Striping patterns are visible, especially in the T_obs uncertainty map; these are likely to be coming from the IRAS data. The uncertainty of T_obs is of the order of 1–3% in bright areas, with a noticeable increase in the inner Galaxy and rising to 5–8% in the most diffuse areas of the sky. The same general trend is seen for τ₃₅₃ but with higher values: 2–5% in bright areas and up to 10% in diffuse areas. The uncertainty of β_obs, based on analysis at 30′ resolution, has a slightly different spatial structure. It is typically of 3–4% with a smaller decrease in bright areas and a noticeable increase in the inner Galaxy to 6–8%.

The reduced χ² of the fit is much smaller than unity over most of the sky, due to the fact that calibration uncertainties are taken into account in the fit to give less weight to data points with less precise calibration (Appendix B). To illustrate this, Fig. 11 shows the distribution function of (Data − Model)/Noise for each frequency used in the fit. The noise used here follows the definition of Eq. (B.1); it takes into account instrumental noise and the uncertainties of the calibration, the zero level, and the CMB subtraction. The range adopted in Fig. 11 corresponds to only ± 1σ. At 353 and 3000 GHz, for most of the sky pixels the data are fitted more tightly (to better than 0.1σ) than at 545 and 857 GHz. This implies that 353 and 3000 GHz have a lot of weight in the estimation of the parameters. The 3000 GHz band provides the only data point on the Wien part of the MBB and therefore strongly influences the determination of T_obs. On the other hand, the 353 GHz band strongly influences the determination of τ₃₅₃ and β_obs because it is the closest to the Rayleigh-Jeans part of the spectrum. The compensating small offsets of the distributions at the other two frequencies might suggest that the adopted model does not adequately describe the data. However, these offsets are well within the calibration uncertainties of the data; the overall reduced χ² is lower than unity and these offsets might be removed by a small systematic change in the relative calibration of the data. Given the actual precision of the calibration, it would be premature to conclude that a more complex model is required to fit the data.

We compared the reduced χ² with that from a fit of the data at 5′ with τ₃₅₃, T_obs, and β_obs as free parameters. We were looking for pixels on the sky for which the two-step fitting procedure provides a reduced χ² greater than unity (i.e., a relatively bad fit) while fitting the three parameters simultaneously at full resolution would provide a better solution with a lower reduced χ². This occurred for only 0.3% of the pixels. These pixels, possibly dominated by galaxies, are grouped in small-scale structures located at high Galactic latitude and away from bright interstellar areas.

3.4. Dust radiance

In the previous sections we have described the properties of the parameters that define the shape of the dust SED. Now we examine the dust radiance or dust integrated intensity defined as $\mathrm{ℛ}\mathrm{=}{\mathrm{\int }}_{\mathit{\nu }}{\mathit{I}}_{\mathit{\nu }}\hspace{0.17em}\mathrm{d}\mathit{\nu }\mathit{.}$(8)Because the grains are in thermal equilibrium, this also represents the energy absorbed. Here we estimate ℛ at each sky position by integrating the MBB fit: $\mathrm{ℛ}\mathrm{=}{\mathrm{\int }}_{\mathit{\nu }}{\mathit{\tau }}_{\mathrm{353}}\hspace{0.17em}{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{T}}_{\mathrm{obs}}\mathrm{\right)}\hspace{0.17em}{\left(\frac{\mathit{\nu }}{\mathrm{353}}\right)}^{{\mathit{\beta }}_{\mathrm{obs}}}\mathrm{d}\mathit{\nu }\mathit{.}$(9)This can be done analytically in terms of the Gamma (Γ) and Riemann zeta functions (ζ): $\mathrm{ℛ}\mathrm{=}{\mathit{\tau }}_{\mathrm{353}}\hspace{0.17em}\frac{{\mathit{\sigma }}_{\mathrm{S}}}{\mathit{\pi }}\hspace{0.17em}{\mathit{T}}_{\mathrm{obs}}^{\mathrm{4}}\hspace{0.17em}{\left(\frac{\mathit{k}{\mathit{T}}_{\mathrm{obs}}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\right)}^{{\mathit{\beta }}_{\mathrm{obs}}}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{4}\mathrm{+}{\mathit{\beta }}_{\mathrm{obs}}\mathrm{)}\hspace{0.17em}\mathit{\zeta }\mathrm{(}\mathrm{4}\mathrm{+}{\mathit{\beta }}_{\mathrm{obs}}\mathrm{)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{4}\mathrm{\right)}\hspace{0.17em}\mathit{\zeta }\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\hspace{0.17em}\mathit{,}$(10)where σ_S is the Stefan-Boltzmann constant, k is the Boltzmann constant, h is the Planck constant, and ν₀ = 3.53 × 10¹¹ Hz. Using the fit parameters described above we produced the all-sky map of ℛ shown in Fig. 12, expressed in units of W m^-2 sr^-1.

Note that even though the calculation of ℛ uses the dust parameters (T_obs, β_obs, τ₃₅₃), ℛ does not suffer from any degeneracy in the fit parameters. In this context the MBB should be seen as an interpolating function; ℛ is not very sensitive to the assumptions made in fitting the SED as long as the fit accounts for the data, including the high-frequency turnover. The uncertainty of ℛ arises mostly from the calibration uncertainty of the data and, to a lesser extent, from the limited number of bands used in the fit⁷.

In thermal equilibrium, ℛ is equal to the amount of light absorbed by dust. Assuming constant properties along the line of sight, including the dust-to-gas ratio, $\mathrm{ℛ}\mathrm{\propto }\mathit{U}\hspace{0.17em}\overline{{\mathit{\sigma }}_{\mathrm{a}}}\hspace{0.17em}{\mathit{N}}_{\mathrm{H}}\mathit{,}$(11)where $\overline{{\mathit{\sigma }}_{\mathrm{a}}}$ is the absorption opacity defined similarly to the emission opacity in Eq. (3), averaged over the size distribution and also, in this case, over the spectrum of the ISRF.

This is complementary to τ₃₅₃, which is also used as a proxy for N_H: ${\mathit{\tau }}_{\mathrm{353}}\mathrm{=}\frac{{\mathit{I}}_{\mathrm{353}}}{{\mathit{B}}_{\mathrm{353}}\mathrm{\left(}{\mathit{T}}_{\mathrm{obs}}\mathrm{\right)}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{e}\hspace{0.17em}\mathrm{353}}\hspace{0.17em}{\mathit{N}}_{\mathrm{H}}\mathit{.}$(12)Division by the Planck function factors out any effects due to spatial variations of the dust temperature (potentially linked to spatial variations of U), but τ₃₅₃ is only proportional to N_H if the dust opacity σ_{e 353} is constant. This limitation does not apply to ℛ, which is independent of σ_{e 353} because of thermal equilibrium; ℛ is simply the energy emitted by dust (Eq. (8)), whatever the shape of the SED and regardless of how efficient the grain cooling is. Thus ℛ is closer to a measured quantity, while τ₃₅₃ is a parameter deduced from a model.

At high Galactic latitudes, where the spatial variations of U and σ_{e 353} are expected to be minimal so that both τ₃₅₃ and ℛ should be proportional to dust column density, comparison of maps of τ₃₅₃ and ℛ reveals another fundamental difference, as illustrated in Fig. 13 for one of the faintest areas in the sky: the map of τ₃₅₃ shows surprisingly strong small-scale fluctuations that are absent in the map of ℛ.

This significant difference is due to the impact of the CIBA, especially its decorrelation in frequency. On the one hand, τ₃₅₃ is the division of I₃₅₃ by B₃₅₃(T_obs) (Eq. (12)) and so is contaminated by the CIBA at not only 353 GHz but also 3000 GHz; i.e., because the 3000 GHz band is the only one in the Wien range, it has a strong weight in the determination of T_obs. Furthermore, the CIBA at 3000 GHz and the CIBA in the Planck bands are weakly correlated, so that T_obs contains most of the information on the CIBA at 3000 GHz. Therefore, through I₃₅₃ and B₃₅₃(T_obs), the map of τ₃₅₃ is affected by the CIBA on both the Rayleigh-Jeans and Wien sides, respectively, resulting in strong small scale fluctuations. On the other hand, because ℛ is obtained by integrating I_ν over frequency, it benefits from the fact that the CIBA decorrelates in frequency; i.e., the integral over frequency of the CIBA is close to zero.

In order to relate Planck dust emission to Galactic reddening (Sect. 6), we also made a fit of the dust model to a version of the Planck and IRAS data from which point sources had been removed (Appendix D). From this fit we have also made maps of τ₃₅₃ and ℛ.

4. The Galactic dust emission observed by Planck

4.1. Spatial variations of the dust parameters and ℛ

The all-sky maps of T_obs, β_obs, and τ₃₅₃ in Fig. 9 – represented as well in a polar orthographic projection in in Figs. 14 and 15 to show details in the high-latitude sky – represent the first attempt to fit these three parameters at the same time over the whole sky. Together, these maps of the fit parameters provide information on the dust SED and, quite likely, on the dust properties and their variations with interstellar environment. They are complemented by the map of ℛ in Figs. 12 and 14. Here we discuss only some broad features of these maps, leaving more detailed analysis to future work.

The mean and standard deviation of T_obs, β_obs, τ₃₅₃, and ℛ are given in Table 3 for several different masks ranked in order of decreasing dust contamination, using σ(τ₃₅₃) as a proxy, and so (mostly) of decreasing sky coverage. These include some masks used in Planck cosmology papers (e.g., Planck Collaboration XV 2014). Note how the ranking is reflected in the means and standard deviations listed.

Over the whole sky, the mean of β_obs is 1.62 and its standard deviation is 0.10. The mean of T_obs is 19.7 K and its standard deviation is 1.4 K. The distribution function of T_obs is slightly positively skewed with a high tail that extends up to 60 K. Only about 100 out of the more than 50 million pixels of the N_side = 2048 map have T_obs< 13 K⁸.

The maps of T_obs and τ₃₅₃ presented here should be compared with the ones published as Planck early results by Planck Collaboration XIX (2011). Apart from the facts that we use a more recent release of Planck data (with a different calibration of the 545 and 857 GHz and with ZE removed) and a slightly different approach to the offset determination⁹, the main difference is that we fit for β_obs while Planck Collaboration XIX (2011) used a fixed value, β_obs = 1.8, a convention shared among all the Planck Early Papers dedicated to dust emission (Planck Collaboration XXV 2011; Planck Collaboration XXIV 2011; Planck Collaboration XXI 2011). Even with these differences in data and methodology, the maps of T_obs are remarkably similar. The map of T_obs presented here is higher by about 1 K than that of Planck Collaboration XIX (2011), due principally to the modification of the calibration of the 545 and 857 GHz channels.

Like in Planck Collaboration XIX (2011), the lowest T_obs values are found in the outer Galaxy and in molecular clouds. In general the well-known molecular clouds have a lower T_obs (15–17 K) and higher β_obs (around 1.8) than in the diffuse ISM. This trend is compatible with the result of Planck Collaboration XXV (2011) who reported a steepening of the SED from diffuse to molecular areas in the Taurus molecular cloud.

Small-scale regions of higher T_obs are seen along the Galactic plane and in many of the Gould Belt clouds, most probably related to the local production of dust-heating photons in Galactic star forming regions. The Magellanic Clouds are clearly visible in the parameter maps with a higher T_obs and lower β_obs (Planck Collaboration XVII 2011).

The main noticeable difference with respect to the early results of Planck Collaboration XIX (2011) is the lower T_obs found here in the inner Galactic plane. This is due to the fact that we fit for β_obs, which appears to have a systematically higher value in the inner Galactic plane, in the range 1.8–2.0. The impact of noise and the CIBA is obviously negligible in this bright area of the sky. The higher β_obs found here clearly provides a better representation of the SED, as shown also by Planck Collaboration Int. XIV (2014). This steepening of the dust SED in the inner Galactic plane is also compatible with the analysis of Herschel observations of that region by Paradis et al. (2012).

One striking feature of the T_obs (polar) map is the increase toward both Galactic poles. Selecting the pixels corresponding to the lowest 1% N_{H i}, the mean T_obs is 20.9 K and the mean β_obs is 1.51. This systematic increase of T_obs was also visible in the early all-sky map of Planck Collaboration XIX (2011) that used a constant β_obs, different offsets, a different 3000 GHz map and no ZE removal for the Planck data. In addition, the values we report for the south Galactic pole mask (mean T_obs = 20.5 K, σ(T_obs) = 1.4 K, mean β_obs = 1.59, σ(β_obs) = 0.13) are compatible with the ones reported by Planck Collaboration Int. XVII (2014) using a correlation method that is insensitive to offsets and ZE removal. The nature of this increase of T_obs over a large scale in the most diffuse areas at high Galactic latitudes, incidentally correlated with lower values of β_obs (see Fig. 15), is still to be understood (see Sect. 5.2) but it is unlikely to be caused by a bias by instrumental noise or the CIBA, which both create small-scale fluctuations.

4.2. β_obs – T_obs relation

The all-sky maps of the fit parameters (Fig. 9) reveal some spatial correlation between the parameters. This is especially clear between β_obs and T_obs as illustrated in Fig. 16, lower, using results for all pixels on the sky. Because it includes so many different regions, this two-dimensional histogram can reveal only global trends, here the general anticorrelation.

This anticorrelation is visible in the faintest parts of the sky, at both small and large scales. It is also seen at the scale of clouds; the Gould Belt clouds have a low T_obs (15–16 K) and high β_obs (~1.8). Several other studies have highlighted similar β_obs − T_obs anticorrelations from observations of specific regions on the sky (Dupac et al. 2003; Désert et al. 2008; Paradis et al. 2010; Planck Collaboration XXV 2011). On the other hand, this behaviour does not extend to the Galactic plane where the two parameters seem to be more correlated than anticorrelated.

As pointed out in Sect. 3.2 (see also Shetty et al. 2009b), instrumental noise is an obvious candidate that might create a β − T anticorrelation. However, the fractional variations of β_obs and T_obs observed here over most of the sky significantly exceed the statistical uncertainties of these parameters taking into account noise and calibration uncertainties (see Fig. 10).

On the other hand, as shown in Sect. 3.2.2 and Appendix C, for faint dust emission the CIBA can produce significant variations of T_obs and β_obs at small scales, and although this effect is in fact observed, it is not accounted for in the error budget. To be quantitative, in the pixels corresponding to the lowest 1% values of N_{H i} the observed standard deviations of these parameters are the largest – σ(T_obs) = 1.7 K and σ(β_obs) = 0.18 – while over the whole sky σ(T_obs) = 1.4 K and σ(β_obs) = 0.10 (see Table 3). Based on the Monte-Carlo simulations presented in Sect. 3.2.3, for values of τ₃₅₃ typical of the faintest 1% pixels of the sky the noise and CIBA produce fluctuations of T_obs and β_obs of the order of δT_obs = 1.7 K and δβ_obs = 0.15, providing a credible explanation for the magnitude of the small-scale variations observed in that mask (Table 3).

Even though noise and the CIBA seem to be responsible for the β − T anticorrelation in the most diffuse areas of the sky, they can cause only small-scale fluctuations because of their flat power spectra. Because the monopole of the CIB was removed from the map, the CIBA does not bias β_obs and T_obs globally on the sky, and cannot produce large-scale variations like the increase of T_obs toward the Galactic poles. We have also checked that these results are largely unaffected by the ZE removal (Appendix A.2).

In brighter regions, our Monte-Carlo simulations (Sect. 3.2.3) indicate that noise and the CIBA introduce variations in T_obs and β_obs (Fig. 8) that are below the observed dispersions (Fig. 16, lower). This is true for more than 90% of the sky. One can appreciate these results by looking directly at the parameter maps (Fig. 9). Away from the most diffuse areas of the sky, where T_obs and β_obs vary at small scale mostly due to the CIBA, the main clouds and interstellar structures that are seen in I_ν and in τ₃₅₃ can be recognized in the maps of T_obs and β_obs.

The broad spectral coverage, the high signal-to-noise (and high signal-to-CIBA) ratio of the data on more than 90% of the sky, and the methodology used to minimize the effects of noise and the CIBA on the fit parameters, combine to produce values of T_obs and β_obs with uncertainties of a few percent, much smaller than their dispersions over the sky (Table 3). We conclude that on most of the sky, the relation between β_obs and T_obs is not an artifact of the data processing (zero levels, ZE correction) or due to noise or the CIBA. This conclusion also holds for the large-scale variations of β_obs and T_obs at high Galactic latitude. On most of the sky, the systematic variations of T_obs and β_obs are related to real changes in the shape of SED of the interstellar dust emission.

Even with data-related effects mitigated, the interpretation of the relationship between the MBB parameters is complex. In particular at this point it is difficult to be definitive about the origin of the relationship between T_obs and β_obs. It depends on details of radiative transfer, of variations in U along the line of sight, and of variations in grain structure and size distribution. To identify the relative roles of dust evolution and line-of-sight integration effects in this observed phenomenon, detailed studies of specific spatially-coherent objects in various interstellar environments and at all scales are needed.

4.3. Dust SED in the diffuse ISM

As described in Planck Collaboration VIII (2014), the calibration scheme for the 545 and 857 GHz data has changed since the Planck Early Results. These channels are no longer calibrated using the FIRAS data, but instead rely on observations of planets as for IRAS, DIRBE, and Herschel. Compared to the previous situation, the calibration factor has been divided by 1.15 at 545 GHz and 1.07 at 857 GHz (Planck Collaboration VIII 2014), so that the specific intensities are now lower.

There are two main impacts on dust modelling. First, the shape of the dust SED is modified, changing the average T_obs and β_obs. The FIRAS average dust SED of the diffuse ISM mask (| b | > 15° and N_{H i}< 5.5 × 10²⁰cm^-2, following the definition of Compiègne et al. 2011) was modelled with T_obs = 17.9 K and β_obs = 1.84 by Planck Collaboration XXIV (2011), compatible with the average SED that they found in selected high Galactic latitude fields using IRAS and the early Planck data¹⁰. With the new calibration, the mean values found for the same mask are significantly different: ⟨ T_obs ⟩ = 20.3 K and ⟨ β_obs ⟩ = 1.59 (see Table 3). The dust parameters found here are similar to those found in external galaxies with Herschel¹¹, even though the Herschel frequency coverage is not as extensive (e.g., Dale et al. 2012).

The second impact is on the value of the dust opacity σ_{e ν} = τ_ν/N_H. The increase in T_obs due to the recalibration lowers τ_ν and the opacity. At 250 μm (1 200 GHz), a reference wavelength often used, Boulanger et al. (1996) obtained σ_{e 1200} = 1.0 × 10^-25 cm² while here for the | b | > 15° mask we obtain σ_{e 1200} = 0.49 × 10^-25 cm² (from σ_{e 353} in Table 4 and β_obs = 1.59).

Changing β_obs directly affects the assessment of the material needed to explain the observed thermal emission; in a MBB fit to the SED, a lower β_obs leads to a higher T_obs and therefore to a lower optical depth, which in turn could be interpreted as a lower column density (or mass), or a lower opacity. We also note that the mean value of β_obs found is lower than used for some components in dust models, like graphite in Draine & Li (2007) where β = 2; when fitting with such a model, a higher radiation field strength U would be needed.

5. Dust emission in relation to gas column density

In the previous section we have described the properties of the parameters that define the dust SED. Now we concentrate on the link between the dust emission and the interstellar gas column density, following on many detailed studies in environments from the diffuse ISM (Boulanger & Pérault 1988; Boulanger et al. 1996) to molecular clouds (Pineda et al. 2008; Goodman et al. 2009).

Here the estimate of gas column density, N_H, accounts for atomic and molecular gas: ${\mathit{N}}_{\mathrm{H}}\mathrm{=}{\mathit{N}}_{\mathrm{H}\hspace{0.17em}{i}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{X}}_{\mathrm{CO}}\hspace{0.17em}{\mathit{W}}_{\mathrm{CO}}\hspace{0.17em}\mathit{,}$(13)where the N_{H i} is from the LAB data assuming optically-thin emission, W_CO is from the Planck¹²CO J = 1 → 0 map (type 3) (Planck Collaboration XIII 2014), and X_CO is not constant but is typically 2 × 10²⁰ H₂ cm^-2 K^-1 km^-1 s (Bolatto et al. 2013). “Dark” neutral matter (Planck Collaboration XIX 2011) is by definition left out in this formulation, though it is among the total that can be traced by γ-rays (Grenier et al. 2005). Ionized gas is left out for lack of a proper template.

5.1. Opacity and dust specific luminosity

The optical depth (τ₃₅₃ here) is often taken as a tracer of N_H but this is only accurate if the opacity is constant (Eq. (3)). This requirement can be assessed in the all-sky map of the opacity σ_{e 353} = τ₃₅₃/N_H in Figs. 17 and 18, smoothed to 30′. Although the large dynamic range over the τ₃₅₃ sky is greatly compressed, so that a linear scale can be used, it is clear that there are changes in opacity, even in the diffuse atomic ISM in the high-latitude sky where N_H is well measured. Related to these changes in opacity are changes in the equilibrium dust temperature (Planck Collaboration XXIV 2011; Planck Collaboration Int. XVII 2014), driving complementary changes in the SED parameter τ₃₅₃ through Eq. (6). This demonstrates how τ₃₅₃ is compromised as a tracer of column density.

We saw in Sect. 3.4 how ℛ compensates for such effects, being smoother than τ₃₅₃. This is expected to carry over into the dust specific luminosity ${\mathit{L}}_{\mathrm{H}}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{ℛ}\mathit{/}{\mathit{N}}_{\mathrm{H}}\hspace{0.17em}\mathit{,}$(14)also shown as an all-sky map in Figs. 17 and 18. At high latitudes this is indeed more uniform. This uniformity and the excursions to both higher and lower values at higher column densities relating to the ambient ISRF are taken up in Sect. 5.2.

Complementing the above, for low-column-density lines of sight with 1 × 10²⁰<N_H< 2.5 × 10²⁰ cm^-2, the dependence of σ_{e 353} on T_obs, and by contrast the relative lack of dependence of L_H on T_obs, are evident in Fig. 19, lower.

The statistics of σ_{e 353} and L_H for the various masks are presented in Table 4 using X_CO = 2 × 10²⁰ H₂ cm^-2 K^-1 km^-1 s where relevant¹². Note how for these normalized quantities the systematic ranking seen in Table 3 is not preserved.

To quantify the trends with column density, Fig. 20 shows the mean and standard deviation of σ_{e 353} and L_H in bins of N_H equally spaced in log. The results for X_CO = [ 1,2,3 ] × 10²⁰ H₂ cm^-2 K^-1 km^-1 s are shown.

5.2. The strength of the ISRF

In thermal equilibrium, ℛ is equal to the amount of light absorbed by dust (see Eq. (11)). In normalized form, ${\mathit{L}}_{\mathrm{H}}\mathrm{\propto }\mathit{U}\hspace{0.17em}\overline{{\mathit{\sigma }}_{\mathrm{a}}}$. Also from Eq. (10) for ℛ evaluated from emission, L_H depends on T_obs, β_obs, and σ_{e 353}. Therefore, under the constraint of thermal equilibrium, the measured ℛ and L_H provide insight into not only U and $\overline{{\mathit{\sigma }}_{\mathrm{a}}}$ relating to absorption, but also their relationship to the SED parameters for emission.

Under the hypothesis of a constant dust-to-gas ratio, constant dust absorption cross section, and constant shape of the ISRF spectrum, i.e., constant $\overline{{\mathit{\sigma }}_{\mathrm{a}}}$, the all-sky map of L_H provides a way to trace the spatial structure of the radiation field U over the whole sky quite directly. The large-scale structure of this map is similar to the map of C⁺/N_{H i}, obtained from lower-resolution data by Bennett et al. (1994), that also traces U.

5.2.1. High latitudes

At high latitudes, best seen in the polar maps, L_H is fairly uniform, much more so than the opacity, as quantified by the relative fractional size of their standard deviations (Table 4). This can also be seen over the low column-density range of Fig. 20, lower, where L_H is constant up to N_{H i} = 5.5 × 10²⁰cm^-2 which is a threshold criterion in the | b | > 15° mask. Even in that mask the standard deviation of L_H is less than 20%. We also note again that this column density is below that for which significant H₂ is seen in the diffuse ISM (Gillmon et al. 2006; Wakker 2006; Rachford et al. 2002, 2009). Furthermore, there is unlikely to be local attenuation of the ISRF at such low column densities¹³. Short of a conspiracy among the several factors affecting L_H, this suggests that each of the factors is fairly uniform in the diffuse atomic high-latitude ISM, up to column densities of at least 5 × 10²⁰cm^-2.

There is a relatively flat trend of L_H with respect to T_obs in Fig. 19, lower, which is for moderate column densities 1 × 10²⁰<N_{H i}< 2.5 × 10²⁰ cm^-2. This uniformity is in contrast to that for σ_{e 353} in the same figure which is anticorrelated with T_obs along a locus of constant L_H, a phenomenon also reported by Planck Collaboration XXIV (2011) and Planck Collaboration Int. XVII (2014). There is also a striking difference between the polar maps of L_H and of T_obs. This demonstrates that T_obs is not a simple tracer of U as is often assumed. In particular, it suggests that the increase of T_obs observed toward the Galactic pole is not a direct result of an increase of U. One interpretation, put forward by Martin et al. (2012), is that grains in different regions of the diffuse ISM retain the effects of different past histories of evolution, e.g., through aggregation and fragmentation, even though the density and timescale argue against present in situ evolution by such processes (Planck Collaboration XXIV 2011). Alternatively, Planck Collaboration Int. XVII (2014) review arguments that grain evolution could be occurring in situ due to UV radiative processing or exposure to cosmic rays. In either case, T_obs would be a response to and tracer of variations in dust properties (grain structure, size distribution, material changes) rather than variations in the strength of the ISRF¹⁴. A corollary is that ℛ could be a better alternative to τ₃₅₃ as a tracer of N_H, at least at high latitudes.

There is a remarkable region near the south Galactic pole with abnormally low L_H, but it does not show up in the ℛ map. Planck Collaboration Int. XVII (2014) argue that it arises because of gas in the Magellanic Stream (MS) that has Galactic velocities and so is counted in N_H while at the same time the dust abundance and dust emission in the low-metalicity MS is very low. Planck Collaboration XXIV (2011) have shown that high velocity clouds (HVC) have relatively low emissivities, which suggests more generally that anomalously low L_H is an interesting diagnostic of HVC-like material that does not have a distinctive HVC velocity. But it is not an argument against using ℛ as a tracer of column density.

However, comparison of the polar maps of IVC (Fig. 2, lower) and L_H shows a correlation of IVC column density with slightly lower L_H over widespread regions. Our interpretation follows Planck Collaboration XXIV (2011) who studied the emissivity of LVC and IVC gas separately and concluded that IVC is Galactic gas that often, though not always, has a lower L_H because dust has been partially destroyed; the dust-to-gas ratio is lower in that component of gas along the line of sight. The amount of such IVC gas would be underpredicted by ℛ¹⁵.

5.2.2. Intermediate to low latitudes

The all-sky map also shows that L_H is not constant, thus strong evidence against using ℛ everywhere as a tracer of column density. For example, the increase of L_H in the inner Galaxy instead implies an increase of the radiation field strength there, by a factor about three compared to the local ISM. The all-sky L_H map also suggests that the ISRF is generally weaker in the outer Galaxy, as expected. Note also that an increase (decrease) in L_H can also be the result of an increase (decrease) in the dust-to-gas ratio accompanying a higher (lower) metallicity.

More localized regions of high L_H are present too: active star formation sites like Cyg X where local sources significantly enhance the IRSF illuminating the dust.

There are localized decreases in L_H as well, coincident with recognizable intermediate-latitude molecular clouds. Our interpretation is that this is a result of a lower ISRF because of attenuation, lowering the energy absorbed by dust within the clouds and hence available to be emitted. In these regions too, ℛ would be compromised as a quantitative linear tracer of column density (see discussion below relating to the Taurus cloud in Fig. 23).

All of these factors contribute to the complicated change of the mean and standard deviation of L_H in Fig. 20, lower, at N_{H i}> 5.5 × 10²⁰cm^-2, the part of the sky that is the complement to the | b | > 15° mask.

5.3. Dust opacity from the diffuse ISM to molecular clouds

The maps of σ_{e 353} in the upper panels of Figs. 17 and 18 reveal variations of the opacity over the sky, variations that are spatially coherent. In the polar plots, the anticorrelation of σ_{e 353} with T_obs is apparent, the same as summarized in Fig. 19 (upper). Because of the β_obs − T_obs anticorrelation discussed in Sect. 4.2, also seen clearly in the polar maps, there is a correlation of σ_{e 353} and β_obs as well. Thus at high latitude τ₃₅₃ is not a reliable measure of N_{H i}.

The general increase of σ_{e 353} by almost a factor of three toward higher column density can be followed on the all-sky map down to intermediate latitudes, to known molecular clouds (e.g., Taurus, Orion, ρ Ophiuchi). These tend to have lower T_obs and higher β_obs¹⁶.

Figure 20, upper, shows the dependence of σ_{e 353} on N_{H i}. There is a small range at low N_H over which σ_{e 353} is at a minimum and roughly constant. But as shown in Fig. 21, where the slope of a fit of τ₃₅₃ vs. N_{H i} over the same range of N_{H i} corresponds to the same σ_{e 353}, the non-linear increase of τ₃₅₃ (and σ_{e 353}) with N_{H i} sets in at a rather low column density. This in a range where as discussed above L_H is constant and there is no significant molecular hydrogen or H i self-absorption. Thus the increase of σ_{e 353} is real and not a reflex of unaccounted dark gas. Again, this compromises τ₃₅₃ as a measure of N_{H i} in the diffuse ISM.

The opacity continues to increase over the range 3 × 10²⁰<N_H< 1 × 10²¹ cm^-2 reaching a plateau thereafter, with a dependency on the choice of X_CO since the gas is predominantly molecular there. The choice of X_CO = 2.0 × 10²⁰ H₂ cm^-2 K^-1 km^-1 s recommended by Bolatto et al. (2013), results in a flat plateau at about twice the value in the diffuse ISM (dotted line).

5.4. Dust at the lowest column densities

At the lowest column densities (N_H< 1 × 10²⁰ cm^-2) we note an increase of σ_{e 353} and L_H (Fig. 20). This effect is also seen directly in the correlation of I_ν vs. N_{H i} in Fig. 5 where all 857 and 3000 GHz data points fall above the correlation for N_{H i}< 1.0 × 10²⁰ cm^-2. It is also the case at 545 and 353 GHz and it thus propagates into the map of τ₃₅₃ and ℛ. We checked that this effect is independent of the removal of the zodiacal emission. It is also present in E(B − V) /N_H using the E(B − V) map of Schlegel et al. (1998) which is based on DIRBE.

Using correlation studies, Planck Collaboration XXIV (2011) showed that H i is a reliable tracer of dust up to at least N_HI = 2 × 10²⁰ cm^-2, or as discussed in Sect. 5.2.1ℛ is a good tracer of N_H to somewhat higher column densities. This suggests that the excess opacity at the lowest N_{H i} seen here in this pixel by pixel analysis is the signature of dust associated with the warm ionized medium (WIM), i.e., interstellar dust that is mixed with ionized hydrogen, H+, that is not traced by H i emission.

Assuming that dust in the WIM has a similar L_H as in the H i, the WIM gas column density needed to explain the rise at low N_{H i} is N_{H +} ~ 3 × 10¹⁹ cm^-217. However, a constant value of N_{H +} cannot explain the shape of the rise of L_H. The rise is more compatible with N_{H +} + N_{H i} ≈ 1.1 × 10²⁰ cm^-2 suggestive of an increase of the ionization fraction of WNM toward lower N_{H i}. The apparent extra dust emission seen here would then come from diffuse regions where H i is partly ionized.

5.5. Discussion

The comparison of dust emission and gas column density reveals an increase of dust opacity of a factor about two from the diffuse ISM to molecular clouds (Fig. 20, upper). In the translucent transition region (3 × 10²⁰<N_H< 2 × 10²¹ cm^-2) H₂ might start rising in importance before CO and self-absorption could begin to affect the 21 cm line emission. These effects are difficult to quantify from the present data but they cannot be responsible for the systematic difference in σ_{e 353} observed between the H i dominated and CO dominated regimes. The presence of “dark gas”, whether from H₂ not traced by CO or from 21 cm self-absorption, would simply flatten the rising profile of σ_{e 353} vs. N_H, reaching the plateau somewhat later. The fact that at N_H ~ 6 × 10²¹ cm^-2, a column density where CO is thought to be a reliable tracer of N_H, L_H dips to the diffuse ISM value while σ_{e 353} remains at the plateau value (Fig. 20) is also consistent with an increased dust opacity in denser regions.

As discussed further in Sect. 6.2, this increase of σ_{e 353} in denser regions is also seen when N_H is estimated using near-infrared colour excess or star counts (Arce & Goodman 1999; Cambrésy et al. 2001; Stepnik et al. 2003; Planck Collaboration XXV 2011; Martin et al. 2012; Roy et al. 2013). It is generally accompanied by a decrease of T_obs, which is quite challenging to explain just with radiative transfer effects (Ysard et al. 2012). When the gas has become dense, an increased σ_{e 353} might be attributed to an increase of dust emissivity related to dust aggregation/coagulation (Ossenkopf & Henning 1994; Ormel et al. 2011; Köhler et al. 2012)¹⁸. However, we have seen opacity changes in the diffuse high-latitude ISM as well that need alternative interpretation if the evolution is in situ (Sect. 5.2.1).

Because ℛ is less affected by the CIBA and because of its correlation with N_H over a larger range in column density (see L_H in Fig. 20, lower), we conclude that ℛ is preferred over τ₃₅₃ as a tracer of column density in the high-latitude diffuse ISM, at least for N_H< 5 × 10²⁰ cm^-2. However, this preference does not hold in molecular clouds and star forming regions where ℛ traces not only the column density but also variations of the radiation field strength due to attenuation and/or local sources of heating photons. In such regions, τ₃₅₃ is the preferred tracer of column density, to the extent that σ_{e 353} is constant there¹⁹. This is supported empirically by the good correlation with the colour excess E(J − K_s) discussed below in Sect. 6.2. However, finding the absolute, rather than relative, column density depends on proper calibration of the opacity σ_{e ν}, which appears to vary with column density and be larger in these regions. Furthermore, caution is advised because the opacity changes from diffuse to dense regions, which might occur over the range of column densities encountered in the region being analysed.

6. Dust emission in relation to extinction

A quantity often used to estimate interstellar column density is visible or near-infrared extinction measured along lines of sight to point sources: stars, globular clusters, galaxies, or quasars. It has been established long ago that there exists a correlation between gas and dust column densities, in particular through the comparison of 21 cm emission and visible extinction (e.g., Lilley 1955). The linear relationship between N_{H i} and E(B − V) was established in the 1970s (Savage & Jenkins 1972; Knapp & Kerr 1974; Ryter et al. 1975; Bohlin et al. 1978)²⁰. Knapp & Kerr (1974) advocated using 21 cm observations of N_{H i} as a proxy for extinction and this correlation, especially as calibrated in the diffuse ISM using measurements on extragalactic objects, has been key to correct extragalactic observations for Galactic reddening.

We have seen in Sect. 5 how the amount of dust emission is, not surprisingly, also correlated with N_H. However, because dust is the agent in both extinction and emission, we make a direct comparison of these observables rather than using N_H as an intermediary. An important example of this direct approach is the proposal by Schlegel et al. (1998) to use dust optical depth obtained from FIR emission (IRAS and DIRBE), rather than H i, to estimate reddening (E(B − V)_SFD), through a correlation calibrated on reddening measurements of galaxies. Such an approach is pursued in Sect. 6.1. For higher column density lines of sight, we compared dust emission to colour excess measurements based on 2MASS stellar photometry²¹.

6.1. Correlation with E(B – V) from quasars

Here, based on the Planck dust emission, we develop a map of E(B − V) applicable to the diffuse ISM at high Galactic latitude. In the years since the work of Schlegel et al. (1998), many models have been put forward self-consistently describing dust emission and extinction (e.g., Draine & Li 2007; Compiègne et al. 2011) and they could be used to convert emission to extinction. However, to be independent of any assumption about dust properties, we decided to remain with an empirical approach. We estimate the conversion factor to E(B − V) using measurements of extinction of extragalactic objects rather than stars to avoid potential biases due to background dust emission.

In particular we estimated E(B − V) using Sloan Digital Sky Survey (SDSS) measurements of quasars. We used the final edition of the SDSS-II quasar catalogue (Schneider et al. 2010) based on the seventh SDSS data release (Abazajian et al. 2009). The catalogue contains 105 783 objects spread over 8 400 deg² mostly on the northern Galactic hemisphere. For each quasar, the observed magnitudes in bands u, g, r, i, and z are given together with their uncertainties. All objects in this catalogue have highly reliable redshift estimates. We limited the sample to a subset of 53 399 quasars at redshifts for which Lyα does not enter the SDSS filters. One benefit compared to the work of Schlegel et al. (1998) is the much larger number of objects. Another is that many studies based on SDSS data have shown that the shape of the extinction curve in the diffuse ISM is compatible with that for stars from Fitzpatrick (1999) with R_V = 3.1 (Jones et al. 2011; Schlafly & Finkbeiner 2011; Mörtsell 2013), so that we can take advantage of all of the multi-colour measurements. The details of how we estimate E(B − V) for each quasar are given in Appendix E.

The correlations of E(B − V) with Planckℛ and τ₃₅₃, from which point sources have been removed (Sect. 3.4), are shown in Fig. 22²². Each is strongly correlated: E(B − V)/ℛ = (5.40 ± 0.09) × 10⁵ and E(B − V) /τ₃₅₃ = (1.49 ± 0.03) × 10⁴. The fractional uncertainty of the slope of the correlation with ℛ is about 20% lower than that for τ₃₅₃. This is not unexpected, both from our discussion of Fig. 13 and, in Sect. 5.2.1, of the factors that influence τ₃₅₃ but not ℛ. This leads us to prefer the solution based on ℛ, E(B − V)_ℛ, to that based on τ₃₅₃, E(B − V)_τ353, for the low column density regions of the sky (see also Sect. 5.5).

The product of E(B − V)/ℛ and the diffuse ISM estimates of L_H gives the ratio E(B − V) /N_H = (1.42–1.46)× 10^-22 mag cm², a factor just 0.82–0.85 lower than that measured using background stars for lines of sight with considerably larger E(B − V) (Bohlin et al. 1978; Rachford et al. 2009). Given all the potential for differences, this agreement is remarkable. On the other hand, the same ratio found using E(B − V) /τ₃₅₃ and the low value of σ_{e 353} in Figs. 20, upper, and 21 and in Table 4 results in a factor 0.55–0.60 lower. Using the E(B − V) map derived from dust emission τ₃₀₀₀ by Schlegel et al. (1998), Liszt (2014) also found a lower ratio, by a factor 0.7, for lines of sight with E(B − V) < 0.1.

6.2. Comparing to E(J – K_s) from star colours in molecular clouds

We have investigated how to estimate E(B − V) for nearby molecular clouds from Planck dust emission. A point of comparison for such regions is mapping of colour excesses E(J − H) and E(H − K_s), or their sum E(J − K_s), based on stellar colours (e.g., Goodman et al. 2009). Here we use maps produced with the AvMAP technique and based on colour excesses from the 2MASS data base (Schneider et al. 2011). Similar maps can be obtained with the NICER and NICEST techniques (e.g., Lombardi et al. 2011). The effective angular resolution of the 2MASS extinction maps used here is a few minutes of arc (see, e.g., Roy et al. 2013), close to that of Planck.

Although the optical extinction is rarely directly measured in the same high column density regions (Martin et al. 2012; Roy et al. 2013), these near-infrared colour excesses are usually expressed as E(B − V)_2MASS after conversion using an assumed shape of the extinction curve, that for stars with R_v = 3.1. This conversion might be inappropriate, and even variable across a field, because of dust evolution affecting all of the colour excess ratios for a given column of dust. Nevertheless, it is still very interesting to compare the spatial details of the dust column density revealed by dust extinction and by dust emission, because each is affected by different systematic effects.

For the Taurus and ρ Ophiuchi molecular clouds, Figs. 23 and 24 present a comparison of four different estimates of E(B − V): E(B − V)_2MASS, E(B − V)_τ353, E(B − V)_SFD, and E(B − V)_ℛ. The Pearson correlation coefficients of E(B − V)_2MASS with the three other maps are given in Table 5. The correlation of E(B − V)_τ353 map with E(B − V)_2MASS is excellent (Pearson coefficient of 0.86 for Taurus and 0.95 for ρ Ophiuchi); this agreement is remarkable given that very different methods and data sets were used to build these two maps. However, notice how the brightest filamentary structures appear with more contrast in E(B − V)_τ353 than in E(B − V)_2MASS, a point to which we shall return.

The Schlegel et al. (1998) map E(B − V)_SFD was also produced from a dust optical depth map (at 3000 GHz in that case; see Sect. 7.3 for details) and so the scale is similar. However, the correlation coefficient with E(B − V)_2MASS is lower and it is clear that a lot of spatial detail is absent. This arises because they estimated T_obs using the lower resolution DIRBE data, thus missing the dust temperature decrease at small scales that accompanies the increase of column density in molecular clouds, and so their map of optical depth which underlies E(B − V)_SFD has lower contrast as well as lower resolution. Being able to follow the small-scale variations of T_obs appears to be essential to gauge properly the full structural details of the molecular clouds. This was revealed by previous studies (Cambrésy et al. 2001; Stepnik et al. 2003; Planck Collaboration XXV 2011), in particular using higher resolution Herschel data (Battersby et al. 2011; Roy et al. 2013), and is now confirmed and reinforced by our Planck analysis.

Although a lot of spatial detail is faithfully reproduced in E(B − V)_ℛ, thanks to the Planck resolution, the correlation with E(B − V)_2MASS is less good. Furthermore, the scale is off. In these denser regions of the ISM, the radiation field strength, and hence ℛ, varies locally due to attenuation and/or local production of photons. The first effect (attenuation) is apparent in Taurus where the densest and brightest filamentary structures appear with less contrast in E(B − V)_ℛ than in E(B − V)_2MASS, and these regions are, consistently, also colder. The case of ρ Ophiuchi is different; because of active star formation, and thus local sources of heating photons, it is a photon-dominated region. The spatial structure of ℛ is therefore visually different than that of τ₃₅₃ or E(B − V)_2MASS, due to the spatial variation of the radiation field.

Opposite to what was found in the diffuse ISM, τ₃₅₃ appears preferable to ℛ as a tracer of column density, in this case E(B − V). But as discussed next it is a complex situation warranting caution. A high correlation coefficient is an important criterion, but the scale and dynamic range are also important.

These effects can be appreciated by the quantification in Fig. 25 where the ratios E(B − V)_τ353/E(B − V)_2MASS and E(B − V)_ℛ/E(B − V)_2MASS are plotted as a function of T_obs.

In Taurus, E(B − V)_τ353 agrees with E(B − V)_2MASS over most of the map, although it is systematically high by typically 25%. Arce & Goodman (1999) reported a similar result in their study of E(B − V)_SFD in Taurus. The exception here is in the coldest parts of the cloud where E(B − V)_τ353/E(B − V)_2MASS increases even more. The systematic departure and spatial variations of E(B − V)_τ353/E(B − V)_2MASS appear to be the result of an increase in σ_{e 353} even in the relative diffuse parts of the map and even more in the coldest (densest) regions. The opacity changes at higher column densities are argued to be related to dust evolution (Planck Collaboration XXV 2011), and unless independently characterized these opacity changes compromise the interpretation of E(B − V)_τ353 as a quantitative measure of dust column density. This is a general concern for all column densities derived from FIR and submillimetre optical depth.

In Taurus, E(B − V)_ℛ also agrees with E(B − V)_2MASS over large parts of the map, although it is slightly low systematically (Fig. 25). In the densest regions, which are cold because of attenuation of the ISRF, L_H is depressed even further and E(B − V)_ℛ/E(B − V)_2MASS decreases. This greatly reduces the contrast across the map of E(B − V)_ℛ.

In ρ Ophiuchi, there is a similar scale difference in the typical E(B − V)_τ353/E(B − V)_2MASS and an upturn toward lower T_obs. Unlike in Taurus, L_H is not generally depressed. There is a correlation of E(B − V)_ℛ/E(B − V)_2MASS with T_obs as expected if T_obs reflects changes in the strength of radiation. Without an independent measure of changes in the ISRF, ℛ is not a reliable quantitative tracer of the dust column density.

6.3. Discussion

As in the comparison of column density measures from dust emission with gas column density, the comparison with dust extinction leads us to the conclusion that ℛ is a slightly better tracer of E(B − V) for diffuse low column density lines of sight. This tracer, of particular interest for extragalactic studies, is the product that we call E(B − V)_xgal in the Planck Legacy Archive (PLA; Appendix F). In addition to its usefulness in estimating Galactic E(B − V) for extragalactic studies, it can also be used to study the structure of the diffuse ISM in regions where the CIBA is dominating the fluctuations at small scales. We stress that E(B − V)_xgal should not be used for E(B − V) > 0.3, where attenuation effects on ℛ become important. In particular, this counter-indication applies in molecular clouds and star forming regions where ℛ traces not only the dust column density but also changes in the ISRF arising from attenuation and/or local sources of radiation. In line with the discussion in Sect. 5.5, in such regions E(B − V)_xgal should not be used.

However, in such regions τ₃₅₃ is well correlated with E(B − V)_2MASS, suggesting an alternative tracer. Again there is an issue with absolute amount, as in the conversion of τ₃₅₃ to N_H, relating to changes in σ_{e 353}. Adopting the scaling factor estimated using the correlation with quasars for more diffuse lines of sight (E(B − V) /τ₃₅₃ = 1.49 × 10⁴ – see Fig. 22) along with the τ₃₅₃ maps could systematically overestimate E(B − V), and this could be exacerbated in the most dense regions where σ_{e 353} increases further.

7. Planck dust products and comparisons with forerunners

7.1. Description of products

As described more fully in Appendix F, the following maps are available in the PLA: the three MBB parameters T_obs, β_obs, and E(B − V)_xgal, together with the associated uncertainty maps. The map of dust integrated intensity, ℛ, can be obtained readily from the three MBB parameter maps using the analytical expression of Eq. (10).

E(B − V)_xgal, a scaled version of ℛ, was obtained from MBB parameters of a fit to data from which point sources have been removed (see Sect. 3.4 and Appendix D). The map of this ℛ can therefore be obtained by dividing the E(B − V)_xgal map by the conversion factor from Sect. 6.1 (Fig. 22, upper).

The main limitations on the Planck dust products from the data themselves are related to the IRAS data. There is some residual striping in the 3000 GHz IRAS data that propagates mainly into the map of T_obs. There is also about 4% of the sky that was not observed by IRAS (Beichman et al. 1988). This area was filled with DIRBE data (Miville-Deschênes & Lagache 2005). Finally, we stress that the dust model is based on data that unavoidably include the CIBA.

7.2. Modelling dust emission: comparison of Planck with Finkbeiner et al. (1999)

One of the important applications of these parameter maps will be to combine them using Eq. (6) to model the SED of the dust emission in the submillimetre range.

A benchmark for comparison is the parametric SED model developed by Finkbeiner et al. (1999). Motivated by data from the FIRAS experiment, they modelled the dust emission as the sum of two MBBs. Each MBBs component is in principle characterized by three parameters. However, by adopting the view that some parameters (or related ratios) are global and by using the constraints provided by thermal coupling to the same radiation field at a given sky position (like Eqs. (10) and (11)) so that the two temperatures are coupled, their model is simplified to only two degrees of freedom instead of six. Thus over the whole sky, fitting IRAS and lower resolution DIRBE data, the model can be summarized by two templates corresponding to the two degrees of freedom: the total dust optical depth at 100 μm at a resolution of 6.́1 and a dust temperature map (for either component) at a resolution of several degrees (it is almost constant at high Galactic latitude).

The dust model can be improved significantly by exploiting the Planck data. The exploration of the parametrization of the dust SED done previously at 7° resolution with FIRAS data can be done at 5′ with much better sensitivity. With the recalibration of the Planck 545 and 857 GHz data we have shown that a single MBB is a good representation of the dust SED over the 353–3000 GHz frequency range, well within the relatively large calibration uncertainties of the data (about 10% at 545, 857, and 3000 GHz) (Sect. 3.3, Fig. 11). Planck Collaboration Int. XVII (2014) reached the same conclusion. Even though the Planck dust model assumes that the dust emission can be modelled by a single MBB from 353 to 3000 GHz, it has one extra degree of freedom compared to Finkbeiner et al. (1999) because τ₃₅₃, T_obs, and β_obs are estimated at each sky position. Finally, the parameter maps are at higher resolution (5′, 5′, and 30′, respectively) and are less noisy, providing a very tight description of the data over that frequency range. With frequency coverage spanning the SED, we can also measure the radiance ℛ.

Figure 26 shows the ratio of the predicted brightness at 353 GHz from the Finkbeiner et al. (1999) dust model (model 7) and that from the Planck model, both at 30′ resolution.

Although the global ratio of the two maps is compatible with one, there are variations at all scales much larger than the uncertainties of the Planck model. Local variations larger than 30% are seen all over the sky, especially in the Galactic plane; the outer Galaxy is significantly underpredicted in the Finkbeiner et al. (1999) model while the inner Galaxy is too bright. Because the Planck dust model is a particularly tight representation of the Planck 353 GHz data (Fig. 11), the same discrepancies are seen by comparing the Finkbeiner et al. (1999) model, integrated in the Planck bandpass, directly with the 353 GHz Planck data.

The Planck dust model produces an accurate 353 GHz map almost free of instrumental noise (but recall that the model includes the effects of the CIBA). That model map along with T_obs can be the basis for extrapolation to lower frequencies, assuming that the appropriate β_obs,mm can be identified.

7.2.1. Frequency range of application

We recall that our fit was done using data from 353 to 3000 GHz. Extrapolating the model outside this range is not recommended. At higher frequencies the dust emission is known to be in excess with respect to the big grain MBB, due to the emission from smaller, stochastically heated, grains (see, e.g., Draine & Li 2007; Compiègne et al. 2011).

At frequencies below 353 GHz the dust SED seems to be flatter than that found for the frequency range here, i.e., from the tests that we have made, extrapolation of the dust model underpredicts the unmodelled Planck dust emission at lower frequencies. This is in accord with the results of Planck Collaboration Int. XVII (2014) for the south Galactic pole area where β_obs,mm = 1.54 ± 0.03 between 100 and 353 GHz and β_obs,FIR = 1.65 ± 0.10 at higher frequencies. The spectral index of dust between 100 and 353 GHz estimated over the whole sky using the Commander-Ruler code (Planck Collaboration XII 2014) has a mean value of 1.49, thus also significantly flatter, and interesting variations (their Fig. 16) that are similar though not identical to what we see in our map of β_obs (Fig. 9, lower). While further discussion is beyond the scope of this paper, it is clear that extrapolating the current model to frequencies lower than 353 GHz needs to be approached with caution.

7.3. Extinction: comparison of Planck with Schlegel et al. (1998)

One of the expected uses of the Planck dust products presented here is to estimate reddening for extragalactic studies. Here we evaluate how the PlanckE(B − V)_xgal map compares with the widely-used E(B − V)_SFD map from Schlegel et al. (1998).

The E(B − V) map of Schlegel et al. (1998) is proportional to their map of τ₃₀₀₀, obtained assuming a constant β_obs and with T_obs estimated with low resolution, low sensitivity DIRBE data. The proportionality factor E(B − V) /τ₃₀₀₀ was estimated by correlating τ₃₀₀₀ with colour excess measurements on 389 galaxies, assuming R_V = 3.1 and the extinction curve from Cardelli et al. (1989) and O’Donnell (1994). This factor E(B − V) /τ₃₀₀₀ has been checked by different probes of extinction and for much larger samples. These studies all used SDSS data and showed that, globally, for E(B − V) < 0.5, the regime of interest for extragalactic studies, E(B − V)_SFD is precise to 15%, though not a fully consistent picture. Using reddening of quasars, Mörtsell (2013) concluded that E(B − V)_SFD underestimates E(B − V) by 20% at low E(B − V) values. Using reddening measurements of elliptical galaxies Peek & Graves (2010) found local variations of E(B − V)_SFD/E(B − V)_elliptical but no systematic bias of the E(B − V)_SFD map. Using reddening of stars, Schlafly et al. (2010) and Schlafly & Finkbeiner (2011) concluded that E(B − V)_SFD overestimates E(B − V) by 14% but with spatial variations of the normalization of the order of 10% which might be attributed to biases in the dust temperature map of Schlegel et al. (1998). Jones et al. (2011) used SDSS colours of M dwarfs to estimate the Galactic extinction properties. In their comparison to Schlegel et al. (1998) they often find E(B − V) values lower than E(B − V)_SFD at lower Galactic latitudes but mention that this could be due to the fact that extinction toward stars does not trace the full line-of-sight dust column density, an effect that could be present in the analysis of Schlafly et al. (2010); Schlafly & Finkbeiner (2011). At high Galactic latitudes Jones et al. (2011) report a number of lines of sight where E(B − V)_SFD seems to underestimate E(B − V).

Based on near-infrared data in brighter areas (E(B − V) > 1.5), Arce & Goodman (1999) and Cambrésy et al. (2001) found that E(B − V)_SFD overestimates E(B − V) systematically by more than 30%. One possible interpretation is an increase of the dust emission efficiency, σ_a, relative to the dust absorption cross section, $\overline{{\mathit{\sigma }}_{\mathrm{a}}}$, in the dense medium, potentially caused by changes in the grain structure, similarly to what is seen with τ₃₅₃/N_H in Fig. 20, upper. All of these studies reveal variations of the ratio E(B − V)_SFD/E(B − V)_SDSS that depend on position on the sky, on methodology (stars, galaxies, quasars), and on the range of E(B − V). There seem to be systematic trends where E(B − V)_SFD overestimates reddening in dense regions and underestimates reddening in diffuse areas at high Galactic latitudes.

The PlanckE(B − V)_xgal and E(B − V)_SFD are compared in Fig. 27, for the low N_{H i} mask (Fig. 3) corresponding to N_{H i}< 2 × 10²⁰ cm^-2 and low IVC emission. The comparison is done after smoothing E(B − V)_ℛ to 6.́1, the resolution of E(B − V)_SFD. The correlation between the two maps is excellent with E(B − V)_SFD = 0.92 E(B − V)_ℛ − 0.003. The dispersion around the correlation is 7%. In relative terms E(B − V)_SFD underestimates E(B − V)_xgal by 8% but in absolute terms the underestimate is more as there is also a negative offset of 0.003 mag.

Our result depends on our choice of using E(B − V) of quasars to calibrate ℛ. Given the variety of methodologies used in the previous studies that have examined the calibration of the high latitude extinction (Schlafly et al. 2010; Schlafly & Finkbeiner 2011; Peek & Graves 2010; Jones et al. 2011; Mörtsell 2013), each having their own biases, the agreement is excellent.

8. Conclusion

We have presented an all-sky model of dust emission based on Planck and IRAS data at 5′ resolution, covering the frequency range from 353 to 3000 GHz. We fit the data at each pixel of the HEALPixN_side = 2048 grid assuming a modified blackbody model. The parameters were estimated using a χ² minimization, in two steps. First the data were smoothed to 30′ and τ₃₅₃, T_obs, and β_obs were estimated in a three-parameter fit. Then, using the spectral index β_obs at 30′ as a fixed parameter, τ₃₅₃ and T_obs were estimated at 5′. We showed that this method minimizes the effect of noise on the determination of the parameters (especially on the T_obs – β_obs degeneracy in faint parts of the sky).

Over the whole sky, the mean and standard deviation of T_obs and β_obs are [19.7, 1.4] K and [1.62, 0.10], respectively. The uncertainties of each parameter are about 3–6%. We showed using Monte Carlo simulations that the uncertainties are dominated by the CIBAs, which are highly correlated within the Planck bands but only at a 30% level between Planck and IRAS 3000 GHz. The CIBA is the dominant source of uncertainty after smoothing the data to 30′ due to its non-white power spectrum (C_ℓ ∝ l^-1).

This Planck dust model reproduces the Planck data well within the noise level at all frequencies. Comparison of the Finkbeiner et al. (1999) dust model with the new data and model shows only broad agreement, with variations of the order of 30% at all scales.

We found an increase of the dust opacity, τ₃₅₃/N_H, by a factor of two from the diffuse to the higher column density (denser) ISM. Empirically, this increase is associated with a decrease in T_obs; because grains are in equilibrium with the interstellar radiation field, we interpret this as a response to the increased dust emissivity. We also noted an excess of dust emission and opacity at H i column densities lower than 10²⁰ cm^-2 that might be attributed to dust in the warm ionized medium (WIM).

The combination of Planck and IRAS data allowed us to model the dust emission from the Rayleigh-Jeans regime over the peak of the SED to the Wien side and therefore to estimate the dust integrated specific intensity or radiance ℛ for each line of sight. We also presented a map of the specific dust luminosity L_H by normalizing ℛ with respect to H i (Eq. (14)). Given thermal equilibrium emission, this is a direct tracer of U, the average strength of the interstellar radiation field (weighted by dust absorption opacity) on each line of sight. This L_H map reveals an increase of U in the inner Galaxy, in active star forming regions, and in the Magellanic Clouds. The map of L_H was shown to be notably different from that of T_obs, indicating that T_obs is not a simple tracer of U as often assumed. This is especially true at high Galactic latitudes where it was found that L_H is fairly uniform and T_obs depends (inversely) on the opacity τ₃₅₃/N_{H i}, confirming early Planck results (Planck Collaboration XXIV 2011). This reveals that τ₃₅₃ is not the most reliable estimator of column density in the diffuse ISM. The spatial variations of T_obs observed in the high-latitude sky appear to be a response to variations of the dust emission opacity resulting in grains of different equilibrium temperature even when exposed to the same U. The analysis at high Galactic latitude is consistent with U being fairly uniform, so that ℛ is a good estimator of column density and can be used to estimate E(B − V) there.

On the other hand, in molecular clouds we showed that variations of T_obs are dominated by the effect of attenuation of the interstellar radiation and/or local sources of heating photons. In this type of environment, where the amplitude of the CIBA is negligible, τ₃₅₃ is a better estimator of column density than ℛ, but the scale depends on the adopted opacity. Compared to the lower resolution work of Schlegel et al. (1998), the MBB analysis of Planck data in this paper provides estimates of T_obs at 5′ resolution and thus an improved higher-resolution estimate of τ₃₅₃, especially in high-contrast molecular regions where the dust temperature and column density vary markedly at small scales.

The Planck dust model was used to produce a map to measure Galactic dust reddening for extragalactic studies at high Galactic latitude, E(B − V)_xgal. This map was based on the radiance ℛ and calibrated by comparison with SDSS reddening measurements of quasars. The correlation of E(B − V)_xgal with the E(B − V) map of Schlegel et al. (1998) is very tight for N_{H i}< 2 × 10²⁰ cm^-2, but has a slope significantly different than one, in the sense that the E(B − V) map of Schlegel et al. (1998) underestimates E(B − V) by 8% in the diffuse ISM. We stress that E(B − V)_xgal is reserved for extragalactic studies; it should not be used to estimate reddening in lines of sight where E(B − V) > 0.3, i.e., where attenuation effects on ℛ become important. There we recommend the map of τ₃₅₃ multiplied by the E(B − V) /τ₃₅₃ ratio also calibrated using quasars. However, systematic decreases of scale can arise from region to region, and even locally within a region, because of the increases in the opacity σ_{e 353} that, empirically, accompany increase in (column) density.

The Planck dust products (τ₃₅₃, T_obs, β_obs, ℛ and E(B − V)_xgal) are available on the PLA (see Appendix F).

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement No. 267934. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the US Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbuk agakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington. Some of the results in this paper have been derived using the HEALPix package.

References

Appendix A: Zodiacal emission

Appendix A.1: Zodiacal emission correction at 100 μm

Zodiacal emission (ZE) is a component that is difficult to remove from the data as it is changing with sky position as well as time of observation. A detailed model has been built to correct the Planck maps for ZE (Planck Collaboration XIV 2014). For the IRAS 100 μm map, as mentioned in Sect. 2.2, the IRIS and SFD maps were not corrected for ZE in the same way. The impact of the different ZE correction for these two maps can be appreciated from Fig. A.1 where the dust emissivity I_ν/N_{H i} is plotted for the four frequencies, using data smoothed to 1° and selecting only pixels with N_{H i}< 3 × 10²⁰ cm^-2. In such a low column density regime, the correlation between dust emission and H i column density was shown to be tight (Fig. 5). Therefore, any systematic variation of I_ν/N_{H i}, especially at low angular resolution, can reveal emission from components other than H i, for example the warm ionized medium (WIM) or residual ZE. To assess the latter, in Fig. A.1 the emissivity is plotted as a function of ecliptic latitude. For the Planck frequencies the emissivity is almost constant with ecliptic latitude, validating the Planck ZE removal. The SFD data also show a constant emissivity, the exception being the IRIS map that shows a systematic increase toward the ecliptic plane, indicative of residual ZE. This is why we implemented the procedure described in Sect. 2.2.

Appendix A.2: Impact of zodiacal emission correction on dust parameters

To evaluate the impact of the ZE on our analysis, we also compared our results on fitting parameters with those of Planck Collaboration Int. XVII (2014) for the same masked region that they studied, an area of 7500 deg² toward the south Galactic cap. This check is useful because there is a fundamental difference between the two analyses. Here the dust SED is modelled using the observed specific intensity for each pixel independently; results are therefore sensitive to uncertainties in the zero levels of the maps. On the other hand, Planck Collaboration Int. XVII (2014) correlated the dust maps with H i within regions 15° in diameter; they showed that their results are insensitive to the zero level of the maps and to the ZE that is very uniform on 15° scales.

To further evaluate the impact of the ZE correction on our analysis, we have explored different data configurations and compared our results and those of Planck Collaboration Int. XVII (2014) in the south Galactic pole area.

Given the lower-resolution results of Planck Collaboration Int. XVII (2014), we used data smoothed to 60′ on an N_side = 128 grid. The comparison was done using different combinations of maps: IRIS or SFD at 100 μm (at such low resolution, SFD is equivalent to the combined IRIS+SFD map, see Sect. 2.2) and maps with and without ZE removed for Planck. To be compatible with Planck Collaboration Int. XVII (2014), the fit was done using a fixed β_obs = 1.65. The results we obtained on T_obs and τ₃₅₃/N_{H i} for this south Galactic cap region are compiled in Table A.1. The differences in T_obs and τ₃₅₃/N_{H i} between data configurations are limited, within the standard deviation observed over the region. They are especially small for T_obs that shows variations of less than 5% between data configurations. The largest effect is from the removal of the ZE in the Planck data that reduces τ₃₅₃ by 15%. The impact of the choice of IRIS or SFD on τ₃₅₃ is more limited; fitting the data with SFD produces a τ₃₅₃ about 3% higher than with IRIS. Nevertheless and even though there is a general good spatial correlation between the maps of T_obs and τ₃₅₃/N_{H i} obtained with the two methods, care should be taken in comparing them in greater detail. Contrary to the results obtained using a fit of the observed specific intensity, the results of Planck Collaboration Int. XVII (2014) are not sensitive to dust emission associated with the WIM that is not spatially correlated with H i and to the CIBAs. In addition, Planck Collaboration Int. XVII (2014) showed that in the southern Galactic cap area there are H i clouds at local velocities that do not have associated dust emission. These effects produce spatial variations of T_obs and τ₃₅₃/N_{H i} computed with the two methods. Even with these caveats, there is a good agreement between the two analyses. In particular we note that the data configuration combining Planck (ZE removed) together with the SFD map at large scales (the equivalent of the 100 μm map built in Eq. (1)) has ⟨ τ₃₅₃/N_{H i} ⟩, the closest to the values found by Planck Collaboration Int. XVII (2014).

Appendix B: The χ² fit

Adopting a reference frequency ν₀ = 353 GHz, the parameters τ_ν0, T_obs and β_obs (see Eq. (6)) were found at each sky position p by fitting the Planck and IRAS data I_ν(p). The MPFIT χ² minimization routine was used (Markwardt 2009). Colour corrections due to finite bandpasses (Planck and IRAS) were taken into account explicitly at each iteration. To speed up the convergence of the fit, initial estimates of the parameters were provided: β_obs was initialized to 1.65, T_obs was initialized to the colour temperature T_3000−857 obtained from the ratio of I₃₀₀₀/I₈₅₇ and assuming β_obs = 1.65, and τ₃₅₃ was initialized to I₃₀₀₀/ [ B₃₀₀₀(T_3000−857) × (3000/353)^1.65 ] for frequencies in GHz.

The fit takes into account the calibration uncertainty of the data (c_νI_ν(p) – see Table 1), the uncertainty of the CMB removal estimated to be c_νI_CMB(p) where I_CMB is the SMICA map (Planck Collaboration XII 2014), the uncertainty of the offset (δO_ν – see Table 1), and the instrumental noise n_ν(p). For both IRAS and Planck, the instrumental noise is modulated inversely by the square root of the coverage map (number of times a sky pixel p has been observed). All sources of uncertainty are added in quadrature: ${\mathit{N}}_{\mathit{\nu }}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{c}}_{\mathit{\nu }}^{\mathrm{2}}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{c}}_{\mathit{\nu }}^{\mathrm{2}}\hspace{0.17em}{\mathit{I}}_{\mathrm{CMB}}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}\mathit{\delta }{\mathit{O}}_{\mathit{\nu }}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}{\mathit{n}}_{\mathit{\nu }}\mathrm{(}\mathit{p}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathit{.}$(B.1)

Appendix C: Impact of the CIBA on the dust parameters

This section describes the Monte Carlo simulations done to quantify the impact of instrumental noise and the CIBA on the parameters of the MBB fit. Specifically, we simulated a single SED as the sum of dust emission, the CIBA, and instrumental noise, for data smoothed to different angular resolution.

The dust emission was modelled using typical values of dust parameters and taking into account the Planck and IRAS bandpasses. The noise level used at each frequency is the median value found in the low N_{H i} mask (Fig. 3), properly scaled as a function of the angular resolution considered. The CIBA levels in the Planck channels follow the parametrization of Planck Collaboration Int. XVII (2014) (T_CIB = 18.3 K and β_CIB = 1.0) normalized at 857 GHz to the value given by Planck Collaboration XVIII (2011). The CIBA level at 3000 GHz is that of Pénin et al. (2012). The scaling of the CIBA level with resolution was done assuming a power spectrum C_ℓ ∝ ℓ^-1. The noise and CIBA levels for each frequency and each resolution are given in Table C.1.

Because galaxies at different redshifts have their peak emission at different frequencies, the structure of the CIB on the sky is only partly correlated between frequencies. Looking at Planck Collaboration XVIII (2011) and Planck Collaboration XXIV (2011), who showed the residual brightness fluctuations in selected faint patches of the sky after removal of the interstellar dust emission traced by 21 cm emission, one can appreciate visually the level of spatial correlation of the CIBA between frequencies. The fluctuations are strongly correlated in the 857–353 GHz range, slightly less at 3000 GHz. The level of frequency (de)coherence of the CIBA was quantified by Planck Collaboration XXX (2014) for 150 <ℓ< 1000. In the high-frequency channels of Planck they report a strong correlation: 0.95 between 857 and 545 GHz, 0.91 between 857 and 353 GHz, and 0.98 between 545 and 353 GHz. It is weaker between 3000 GHz and the Planck frequencies: 0.36, 0.31, and 0.29 at 857, 545 and 353 GHz respectively. Planck Collaboration Int. XVII (2014) modelled the part of the CIBA correlated with the ones at 857 GHz based on the results of Planck Collaboration XXX (2014). They found that it can be well fitted by a MBB function with T_CIB = 18.3 K and β_CIB = 1.0. This parametrization is compatible with the one found in a similar study by Planck Collaboration XXIV (2011).

To model the CIBA, c_ν, taking into account the inter-frequency correlation, we have proceeded the following way. We made the assumption that c_ν is fully correlated in the Planck frequencies. Therefore, the CIBA in the Planck bands is simply modelled as a single random draw scaled by the CIBA levels: ${\mathit{c}}_{\mathit{\nu }}\mathrm{=}{\mathit{A}}^{\mathrm{cor}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{CIBA}\hspace{0.17em}\mathit{\nu }}\hspace{0.17em}\mathit{,}$(C.1)where A^cor is drawn from a normal distribution with unit standard deviation and zero mean, and σ_{CIBA ν} is the CIBA level at each Planck frequency given in Table C.1.

At 3000 GHz the CIBA is not strongly correlated with the CIBA at Planck frequencies so it was divided into one part correlated with 857 GHz and another part uncorrelated: ${\mathit{c}}_{\mathrm{3000}}\mathrm{=}{\mathit{A}}^{\mathrm{cor}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{CIBA}\hspace{0.17em}\mathrm{3000}}^{\mathrm{cor}}\mathrm{+}{\mathit{A}}^{\mathrm{uncor}}\hspace{0.17em}{\mathit{\sigma }}_{\mathrm{CIBA}\hspace{0.17em}\mathrm{3000}}^{\mathrm{uncor}}\hspace{0.17em}\mathit{,}$(C.2)where A^uncor represent a second random draw. The correlated part, ${\mathit{\sigma }}_{\mathrm{CIBA}}^{\mathrm{cor}}$, is estimated using the [T_CIBA = 18.3 K, β_CIBA = 1.0] parametrization, normalized with c₈₅₇. The uncorrelated part, ${\mathit{\sigma }}_{\mathrm{CIBA}}^{\mathrm{uncor}}$, is simply the quadratic complement to σ_{CIBA 3000}. At 5′ these contributions are ${\mathit{\sigma }}_{\mathrm{CIBA}\hspace{0.17em}\mathrm{3000}}^{\mathrm{cor}}\mathrm{=}\mathrm{0.049}$ MJy sr^-1 and ${\mathit{\sigma }}_{\mathrm{CIBA}\hspace{0.17em}\mathrm{3000}}^{\mathrm{uncor}}\mathrm{=}\mathrm{0.087}$ MJy sr^-1 highlighting the fact that the uncorrelated part is the dominant source of CIB fluctuations at 3000 GHz.

Noise was assumed to be uncorrelated in frequency, with levels given in Table C.1.

Depending on the situation, this procedure was executed with different dust parameters (τ₃₅₃, T_obs, and β_obs) and at different data resolutions. The fit of the simulated SED was done using the method that was used with the real data, including colour corrections. Simulations were done with and without noise to quantify the specific effect of noise and the CIBA on the recovered parameters.

Appendix D: Production of maps with point sources removed

The E(B − V)_xgal product, to be used for estimating Galactic reddening for extragalactic studies, is obtained with a fit of the dust model on a version of the Planck and IRAS data from which point sources have been removed. This section describes how the removal was done. Most of the software used in this process is part of the Planck Sky Model software package (PSM). A description of the PSM and how to download it can be found in Delabrouille et al. (2013).

Appendix D.1: Making the mask containing the radio and IR sources

A source mask is made as a HEALPix map from the union of discs centred on a selected set of known FIR and radio sources. The FIR sources are compiled from the IRAS Point Source Catalog (PSC, Beichman et al. 1988) and the Faint Source Catalog (FSC, Moshir et al. 1992). Radio sources are derived from the WMAP catalog of point sources detected in the 9-year sky maps (Bennett et al. 2013), supplemented with a selection of extragalactic point sources detected by Planck at 100 GHz (ERCSC, Planck Collaboration VII 2011) that are not in the WMAP catalog. Pixels within a 15′ radius of each source are masked.

Appendix D.2: Grouping adjacent pixels for each source

While most sources are isolated, in rare cases masks around two or more sources overlap. This results in larger, non-circular holes in the mask. As a first step, we divide the mask into a set of small connected holes (most but not all of which correspond to the masking of one single source). This is performed by the PSM routine GROUP_ADJACENT_PIXELS.

Appendix D.3: Interpolating across the source

Finally, we estimate the diffuse brightness in each of the small connected holes containing point sources using the PSM routine FILL_SMALLGAP. The interpolation is done using a minimum curvature surface (IDL routine MIN_CURV_SURF). This step produces the final, de-sourced maps.

Appendix E: Estimating E (B – V) from colour excess measurements

This section describes how E(B − V) was estimated from the multi-colour SDSS measurements of quasars. The colour excess due to dust extinction is defined as $\begin{array}{ccc}\mathit{E}\mathrm{(}\mathit{X}\mathrm{-}\mathit{Y}\mathrm{)}& \mathrm{=}& {\mathit{A}}_{\mathit{X}}\mathrm{-}{\mathit{A}}_{\mathit{Y}}\\ & \mathrm{=}& \mathrm{(}{\mathit{m}}_{\mathit{X}}\mathrm{-}{\mathit{m}}_{\mathit{X}\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{m}}_{\mathit{Y}}\mathrm{-}{\mathit{m}}_{\mathit{Y}\mathrm{0}}\mathrm{)}\\ & \mathrm{=}& \mathrm{-}\mathrm{2.5}\log \left(\frac{{\mathit{F}}_{\mathit{X}}}{{\mathit{F}}_{\mathit{X}\mathrm{0}}}\frac{{\mathit{F}}_{\mathit{Y}\mathrm{0}}}{{\mathit{F}}_{\mathit{Y}}}\right)\mathit{,}\end{array}$where, for a given band i, A_i is the extinction, m_i the observed magnitude, m_i0 the absolute magnitude, F_i the observed flux, and F_i0 the intrinsic flux of the source.

The quantity F_i0 corresponds to the source flux density f_λ convolved with the filter transmission T_i(λ)${\mathit{F}}_{\mathit{i}\mathrm{0}}\mathrm{=}{\mathrm{\int }}_{\mathit{i}}{\mathit{T}}_{\mathit{i}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathit{\lambda }}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\hspace{0.17em}\mathit{,}$(E.4)while F_i is the same quantity but affected by extinction; ${\mathit{F}}_{\mathit{i}}\mathrm{=}{\mathrm{\int }}_{\mathit{i}}{\mathit{T}}_{\mathit{i}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathit{\lambda }}\hspace{0.17em}{\mathrm{10}}^{\mathrm{-}\mathrm{0.4}\hspace{0.17em}\mathit{E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{)}\hspace{0.17em}{\mathit{C}}_{\mathit{\lambda }}}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\hspace{0.17em}\mathit{,}$(E.5)where C_λ is the extinction curve for a given value of R_V and normalized to E(B − V): ${\mathit{C}}_{\mathit{\lambda }}\mathrm{=}{\mathit{A}}_{\mathit{\lambda }}\mathit{/}\mathit{E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{)}\hspace{0.17em}\mathit{.}$(E.6)In what follows, we use C_λ of Fitzpatrick (1999) and assume R_v = 3.1.

Equations (E.1) to (E.5) relate the observed colour excess of a source in bands X and Y to its intrinsic spectrum f_λ, to the transmission of the filters T_i and to the extinction of dust along the line of sight. As mentioned by Fitzpatrick & Massa (2005), if the intrinsic spectrum of the source and the filter transmissions are known, the normalization of the extinction curve (i.e., E(B − V)) can be estimated directly from the observed colour excess E(X − Y).

An often overlooked fundamental fact follows from this description. Because the observed magnitudes m_i are obtained with broad band filters, the measured value of E(X − Y) for X = B and Y = V will in general be different from the normalization of the extinction curve E(B − V) in Eq. (E.5). The observed colour excess is indeed a convolution of the source spectrum, the filter transmission and the extinction curve, and so there are bandpass corrections. We stress that the quantity E(B − V) we are estimating here is independent of the spectrum of the background source. It is the value that is used to scale the extinction curve.

From the SDSS quasar catalogue one can deduce directly the observed colours (m_X − m_Y) but to estimate E(B − V) using the formalism described in the previous section one also needs the source spectrum f_λ and the intrinsic colour (m_X0 − m_Y0). One of the advantages of using quasars is that they are known to have a fairly constant spectrum. Below we use the composite quasar spectrum of Vanden Berk et al. (2001) for f_λ. The intrinsic colour (m_X0 − m_Y0) is estimated by correlating (m_X − m_Y) with a tracer of dust extinction. We use the H i column density provided in the SDSS quasar catalogue as a proxy²³: $\mathrm{(}{\mathit{m}}_{\mathit{X}}\mathrm{-}{\mathit{m}}_{\mathit{Y}}\mathrm{)}\mathrm{=}\mathit{\psi }\hspace{0.17em}{\mathit{N}}_{\mathrm{H}\hspace{0.17em}{i}}\mathrm{+}\mathrm{(}{\mathit{m}}_{\mathit{X}\mathrm{0}}\mathrm{-}{\mathit{m}}_{\mathit{Y}\mathrm{0}}\mathrm{)}\hspace{0.17em}\mathit{.}$(E.7)The constant term of this linear regression is the average intrinsic colour of quasars but, as quasars are at different redshifts (up to z ~ 5), this intrinsic colour depends on z. Therefore we correlated (m_X − m_Y) with N_{H i} for quasars in bins of z where the width of each bin is set to have 1000 quasars per bin. Figure E.1 shows the intrinsic colours found as a function of z for the following four colours: (m_u − m_z), (m_u − m_i), (m_g − m_r) and (m_g − m_i).

In order to validate the approach we also computed the intrinsic colour as a function of z of the composite quasar spectrum: $\mathrm{(}{\mathit{m}}_{\mathit{X}\mathrm{0}}\mathrm{-}{\mathit{m}}_{\mathit{Y}\mathrm{0}}\mathrm{)}\mathrm{=}\mathrm{-}\mathrm{2.5}\log \left(\frac{{\mathit{F}}_{\mathit{Y}\mathrm{0}}}{{\mathit{F}}_{\mathit{X}\mathrm{0}}}\right)\mathrm{·}$(E.8)The results for the same four colours are shown in red in Fig. E.1. The good agreement between the two methods indicates that the composite spectrum is a good average representation. It also gives confidence in its use to estimate E(B − V) using the formalism described earlier.

The correspondence between the two methods is less good at high redshift. This can be attributed to the Lyman-α line that enters the shortest wavelength band. Indeed Mörtsell (2013) mentioned that at z> 2.3, Lyman-α introduces large variations in the intrinsic quasar colour computed with band g. A similar effect is seen at z> 1.7 for colours using the u band. We also noticed a larger dispersion in the residual ((m_X − m_Y) − ψN_{H i}) for bins with z< 0.7. Based on this, we selected quasars in the redshift range 0.7 <z< 1.7. The sample contains 53 399 quasars.

For each quasar at a given redshift z, we computed the colour excesses E(X − Y) = (m_X − m_Y) − (m_X0 − m_Y0) using the intrinsic colour estimated with the correlation with N_{H i} in bins of redshift (Eq. (E.7)). We then integrated numerically Eqs. (E.5) and (E.4) using the composite quasar spectrum redshifted appropriately and solved for E(B − V) using Eq. (E.3). Because we work with colours, the exact normalization of the quasar spectrum (f_λ) and of the SDSS transmissions (T_i) cancels out in Eq. (E.3). With this procedure we obtain an estimate of E(B − V) for each quasar independently for each colour excess E(X − Y).

For each quasar position we extracted the values of τ₃₅₃ and ℛ  and then correlated globally with E(B − V) for each colour (X − Y). Table E.1 gives the conversion factors E(B − V) /τ₃₅₃ and E(B − V)/ℛ estimated for each colour as well as the one obtained by combining all colours, weighted by their uncertainties.

Given the very low value of E(B − V), <0.1, the large number of measurements is key here to beat down the noise of the SDSS data but also the variations of the quasar spectra around the template spectrum. The intercepts of the regressions of E(B − V) vs. τ₃₅₃ and E(B − V) vs. ℛ are small for each colour, indicating that there is a coherence between the zero levels of the Planck and IRAS maps used to build τ₃₅₃ and the estimate of the intrinsic colours of quasars, both based on a correlation with N_{H i}.

Appendix F: Maps in the Planck Legacy Archive

Maps of the Planck thermal dust emission model described in this paper can be obtained from the Planck Legacy Archive (PLA)²⁴. The maps available give the three MBB parameters (τ₃₅₃, T_obs, and β_obs) with their uncertainties (δτ₃₅₃, δT_obs, and δβ_obs) obtained using data from which point sources were not removed. In addition, the maps of ℛ (point sources in) and of E(B − V)_xgal are made available; the latter is ℛ obtained with data from which point sources were removed and scaled to E(B − V) with the conversion factor computed using SDSS quasars.

A first release of the Planck thermal dust model was made available on the PLA in March 2013, together with other 2013 Planck products. This first model is superseded by the model presented in this paper. The first model was made using Planck data from which ZE was not removed, to be compatible with what was used for the cosmological analysis. The 3000 GHz map used in this first edition was the IRIS map. The second release described in this paper is based on Planck data from which ZE was removed and on a combination of IRIS and Schlegel et al. (1998) for the 3000 GHz map.

The impact of this change is significant only in regions of the sky of low Galactic emission and low ecliptic latitude. It affects the parameters T_obs and β_obs only slightly, at the 2% level. Considering the pixels in the low N_{H i} mask, the average T_obs went from 20.4 to 20.8 K from the first to the second version, while the average β_obs went from 1.59 to 1.55. That means that the shape of the dust SED is not significantly different from one data set to the other. The main effect of the change of data set is on τ₃₅₃ because the ZE removed data have less emission; from the first to the second version, ⟨ τ₃₅₃ ⟩ is reduced by about 25% in the most diffuse areas of the sky, from 1.3 × 10^-6 to 9.6 × 10^-7.

A significant difference between the two releases concerns the E(B − V) map. In both releases the calibration was done the same way, using SDSS quasars, but in the first release the E(B − V) map was based on τ₃₅₃ instead of ℛ. Globally both E(B − V) maps agree but not on small scales where the 353 GHz the CIBA is much more present in τ₃₅₃ than some spectrally-averaged CIBA is in ℛ. The CIBA in τ₃₅₃ introduces an additional noise in E(B − V) of the order of 0.003 magnitude. In that respect the E(B − V) map of the second version is a significant improvement.

All Tables

All Figures