© ESO, 2016

1. Introduction

Studying the interstellar medium (ISM) is important in a wide range of astronomical disciplines, from star and planet formation to galaxy evolution. Dust changes the appearance of galaxies by absorbing ultraviolet (UV), optical, and infrared (IR) starlight, and emitting mid-IR and far-IR (FIR) radiation. Dust is an important agent in the chemical and thermodynamical evolution of the ISM. Physical models of interstellar dust that have been developed are constrained by such observations. In the present work, we study the ability of a physical dust model to reproduce IR emission and optical extinction observations, using the newly available Planck¹ data.

The Planck data provide a full-sky view of the Milky Way (MW) at submillimetre (submm) wavelengths, with much higher angular resolution than earlier maps made by the Diffuse Infrared Background Experiment (DIRBE; Silverberg et al. 1993) on the Cosmic background explorer (COBE) spacecraft (Boggess et al. 1992). These new constraints on the spectral energy distribution (SED) emission of large dust grains were modelled by Planck Collaboration XI (2014, hereafter Pl using a modified blackbody (MBB) spectral model, parameterized by optical depth and dust temperature. That study, along with previous Planck results, confirmed spatial changes in the dust submm opacity even in the high latitude sky (Planck Collaboration XXIV 2011; Planck Collaboration Int. XVII 2014). The dust temperature, which reflects the thermal equilibrium, is anti-correlated with the FIR opacity. The dust temperature is also affected by the strength of the interstellar radiation field (ISRF) heating the dust. The bolometric emission per H atom is almost constant at high latitude, consistent with a uniform ISRF, but over the full sky, covering lines of sight through the Galaxy, the ISRF certainly changes. The all-sky submm dust optical depth was also calibrated in terms of optical extinction. However, no attempt was made to connect these data with a self-consistent dust model, which is the goal of this complementary paper.

Several authors have modelled the dust absorption and emission in the diffuse ISM, e.g. Draine & Lee (1984), Desert et al. (1990), Dwek (1998), Zubko et al. (2004), Compiègne et al. (2011), Jones et al. (2013), Siebenmorgen et al. (2014). We focus on one of the most widely used dust models presented by Draine & Li (2007, hereafter DL). Earlier, Draine & Lee (1984) studied the optical properties of graphite and silicate dust grains, while Weingartner & Draine (2001) and Li & Draine (2001) developed a carbonaceous-silicate grain model that has been quite successful in reproducing observed interstellar extinction, scattering, and IR emission. DL presented an updated physical dust model, extensively used to model starlight absorption and IR emission. The DL dust model employs a mixture of amorphous silicate grains and carbonaceous grains. The grains are assumed to be heated by a distribution of starlight intensities. The model assumes optical properties of the dust grains and the model SEDs are computed from first principles.

The DL model has been successfully employed to study the ISM in a variety of galaxies. Draine et al. (2007) employed DL to estimate the dust masses, abundances of polycyclic aromatic hydrocarbon (PAH) molecules, and starlight intensities in the Spitzer Infrared Nearby Galaxies Survey – Physics of the Star-Forming ISM and Galaxy Evolution (SINGS, Kennicutt et al. 2003) galaxy sample. This survey observed a sample of 75 nearby (within 30 Mpc of the Galaxy) galaxies, covering the full range in a three-dimensional parameter space of physical properties, with the Spitzer Space Telescope (Werner et al. 2004). The Key Insights on Nearby Galaxies: a FIR Survey with Herschel (KINGFISH, Kennicutt et al. 2011) project, additionally observed a subsample of 61 of the SINGS galaxies with the Herschel Space Observatory (Pilbratt et al. 2010). Aniano et al. (2012) presented a detailed resolved study of two KINGFISH galaxies, NGC 628 and NGC 6946, using the DL model constrained by Spitzer and Herschel photometry. Aniano et al. (in prep.) extended the preceding study to the full KINGFISH sample of galaxies. Draine et al. (2014, hereafter DA14), presented a resolved study of the nearby Andromeda galaxy (M31), where high spatial resolution can be achieved. The DL model proved able to reproduce the observed emission from dust in the KINGFISH galaxies and M31. Ciesla et al. (2014) used the DL model to fit the volume limited, K-band selected sample of galaxies of the Herschel Reference Survey (Boselli et al. 2010), finding it systematically underestimated the 500 μm photometry.

The new Planck all-sky maps, combined with ancillary Infrared Astronomical Satellite (IRAS, Neugebauer et al. 1984) and Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) maps allow us to explore the dust thermal emission from the MW ISM with greater spatial resolution and frequency coverage than ever before. Here we test the compatibility of the DL dust model with these new observations.

We employ WISE 12² (12 μm), IRAS 60 (60 μm), IRAS 100 (100 μm), Planck 857 (350 μm), Planck 545 (550 μm), and Planck 353 (850 μm) maps to constrain the dust emission SED in the range 10 μm <λ< 970 μm. These data allow us to generate reliable maps of the dust emission using a Gaussian point spread function (PSF) with 5′ full width at half maximum (FWHM). Working at lower resolution (1° FWHM), we can add the DIRBE 140 and DIRBE 240 photometric constraints.

We employ the DL dust model to characterize:

• the dust mass surface density Σ_Md;

• the dust optical extinction A_V;

• the dust mass fraction in small PAH grains q_PAH;

• the fraction of the total luminosity radiated by dust that arises from dust heated by intense radiation fields, f_PDR;

• the starlight intensity U_min heating the bulk of the dust.

The estimated dust parameters for M31 are compared with those derived using the independent maps in DA14.

We compare the DL optical extinction estimates with those of Pl-MBB. We further compare the DL model reddening estimates with near-IR reddening estimates from quasi-stellar objects (QSOs) from the Sloan Digital Sky Survey (SDSS, York et al. 2000), and from stellar reddening maps in dark clouds obtained from the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006). These reveal significant discrepancies that call for a revision of the DL model. We find an empirical parameterization that renormalizes the current DL model and provides insight into what is being compensated for through the renormalization.

We use the DL model parameter U_min to bin the Planck and IRAS data for the diffuse ISM and compress the spectral information to a family of 20 dust SEDs, normalized per A_V, which may be used to test and empirically calibrate dust models. We also provide the Planck 217 (1.38 mm) and Planck 143 (2.10 mm) photometric constraints, which are not used in the current dust modelling.

This paper is organized as follows. We describe the data sets in Sect. 2. In Sect. 3, we present the DL dust model, the model parametrization and the method we use to fit the model to the data. The model results are described in Sect. 4. We present the maps of the model parameters (Sect. 4.1), analyse the model ability to fit the data (Sect. 4.2) and assess the robustness of our determination of the dust mass surface density (Sect. 4.3). We compare the dust A_V estimates with the MBB all-sky modelling results from Pl-MBB in Sect. 5. In Sect. 6 we propose an empirical correction to the DL A_V estimates based on the comparison with QSO SDSS data. The DL model is used to compress the Planck and IRAS data into a family of 20 dust SEDs that account for the main variations of the dust emission properties in the diffuse ISM (Sect. 7). In Sect. 8, we extend our assessment of the DL A_V map to molecular clouds. The difference between A_V derived from the DL model and optical observations is related to dust emission properties and their evolution within the ISM in Sect. 9. We conclude in Sect. 10. The paper has four appendices. In Appendix A, we compare our analysis of Planck and IRAS observations of M31 with earlier DL modelling of Herschel and Spitzer data. In Appendix B, we detail how we estimate A_V towards QSOs observed with SDSS. In Appendix C, we analyse the impact of cosmic infrared background (CIB) anisotropies in our dust modelling. In Appendix D, we describe the data products made public in the Planck Legacy Archive.

2. Data sets

We use the full mission maps of the high frequency instrument of Planck that were made public in February 2015 (Planck Collaboration I 2015; Planck Collaboration VIII 2015). The zodiacal light has been estimated and removed from the maps (Planck Collaboration XIV 2014). We remove the cosmic microwave background, as estimated with the SMICA algorithm (Planck Collaboration IX 2015), from each Planck map. Following Pl-MBB we do not remove the CO(3-2) contribution to the Planck 353 GHz band³.

A constant offset (listed in the column marked removed CIB monopole of Table 1) was added to the maps by the Planck team to account for the CIB in extragalactic studies, and we proceed to subtract it. Since the Planck team calibrated the zero-level of the Galactic emission before the zodiacal light was removed, an additional offset correction is necessary. To determine this (small) offset, we proceed exactly as in Sect. 5 of Planck Collaboration VIII (2014) by correlating the Planck 857 GHz map to the Leiden/Argentine/Bonn Survey of Galactic H i, and then cross-correlating each of the lower frequency Planck maps to the 857 GHz map, over the most diffuse areas at high Galactic latitude. These offsets make the intercepts of the linear regressions between the Planck and H i emission equal to zero emission for a zero H i column density. We note that we did not attempt to correct the Planck maps for a potential residual of the CMB dipole as done by Pl-MBB. This is not crucial for our study because we do not use microwave frequencies for fitting the DL model.

We complement the Planck maps with IRAS 60 and IRAS 100 maps. We employ the IRAS 100 map presented in Pl-MBB. It combines the small scale (<30′) features of the map presented by the Improved reprocessing of the IRAS survey (IRIS, Miville-Deschênes & Lagache 2005), and the large scale (>30′) features of the map presented by Schlegel et al. (1998, hereafter SFD). The zodiacal light emission has been estimated and removed from the SFD map, and therefore it is removed from the map we are employing⁴. We employ the IRAS 60 map presented by the IRIS team, with a custom estimation and removal of the zodiacal light⁵. The zero level of the IRAS maps is adjusted so it is consistent with the Planck 857 GHz map.

In Sects. 5 and 9.1, we compare our work with the MBB all-sky modelling results from Pl-MBB. To perform a consistent comparison, in these sections we use the same Planck and IRAS data as in Pl-MBB, i.e. the nominal mission Planck maps corrected for zodiacal light (Planck Collaboration I 2014) with the monopole and dipole corrections estimated by Pl-MBB.

WISE mapped the sky at 3.4, 4.6, 12, and 22 μm. Meisner & Finkbeiner (2014) presented a reprocessing of the entire WISE 12 μm imaging data set, generating a high resolution, full-sky map that is free of compact sources and was cleaned from several contaminating artefacts. The zodiacal light contribution was subtracted from the WISE map assuming that on spatial scales larger than 2° the dust emission at 12 μm is proportional to that in the Planck 857 GHz band. This effectively removes the zodiacal emission but at the cost of losing information on the ratio between WISE and Planck maps on scales larger than 2°. About 18% of the WISE map is contaminated by the Moon or other solar system objects. Aniano (A15, priv. comm.) prepared a new WISE map, with an improved correction of the contaminated area, which we use in this paper.

For typical lines of sight in the diffuse ISM, the dust SED peaks in the λ = 100−160 μm range. DIRBE produced low resolution (FWHM = 42′) all-sky maps at 140 and 240 μm, which can be used to test the robustness of our modelling. Additionally, we perform a lower resolution (1° FWHM) modelling, including the DIRBE 140 and DIRBE 240 photometric constraints. We use the DIRBE zodiacal light-subtracted mission average (ZSMA) maps. This modelling allows us to evaluate the importance of adding photometric constraints near the dust SED peak, which are absent in the Planck and IRAS data.

The FIS instrument (Kawada et al. 2007) on board the AKARI spacecraft (Murakami et al. 2007) observed the sky at four FIR bands in the 50−180 μm range. Unfortunately, AKARI maps were made public (Doi et al. 2015) after this paper was submitted. Moreover, the way their mosaic tiles are chosen and significant mismatch of the zero level of the Galactic emission among the tiles prevent a straightforward integration of the AKARI data into the present modeling.

All maps were convolved to yield a Gaussian PSF, with , slightly broader than all the native resolution of the Planck maps. We use the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) of a sphere coordinates (Górski et al. 2005)⁶. We work at resolution N_side = 2 048, so the maps have a total of 12 × 2 048 × 2 048 = 50 331 648 pixels. Each pixel is a quadrilateral of area 2.94 arcmin² (i.e. about 1.′7 on a side). All maps and results presented in the current paper are performed using this resolution, except those of Sects. 4.3.2 and 7. The most relevant information on the data sets that are used is presented in Table 1. The amplitudes of the CIB anisotropies (CIBA) depend on the angular scale; the values listed in Table 1 are for the 5′ resolution of our data modelling.

3. The DL model

The DL dust model is a physical approach to modelling dust. It assumes that the dust consists of a mixture of amorphous silicate grains and carbonaceous grains heated by a distribution of starlight intensities. We employ the Milky Way grain size distributions (Weingartner & Draine 2001), chosen to reproduce the wavelength dependence of the average interstellar extinction within the solar neighbourhood. The silicate and carbonaceous content of the dust grains has been constrained by observations of the gas phase depletions in the ISM. The carbonaceous grains are assumed to have the properties of PAH molecules or clusters when the number of carbon atoms per grain N_C ≲ 10⁵, but to have the properties of graphite when N_C ≫ 10⁵. DL describes the detailed computation of the model SED, and AD12 describes its use in modelling resolved dust emission.

3.1. Parameterization

The IR emission of the DL dust model is parametrized by six parameters, Σ_Md,q_PAH,U_min,U_max,α, and γ. The definition of these parameters is now reviewed. The model IR emission is proportional to the dust mass surface density Σ_Md. The PAH abundance is measured by the parameter q_PAH, defined to be the fraction of the total grain mass contributed by PAHs containing N_C< 10³ C atoms⁷. As a result of single-photon heating, the tiny PAHs contributing to q_PAH radiate primarily at λ< 30 μm , and this fraction is constrained by the WISE 12 band. Weingartner & Draine (2001) computed different grain size distributions for dust grains in the diffuse ISM of the MW, which are used in DL. The models in this MW3.1 series are all consistent with the average interstellar extinction law⁸, but have different PAH abundances in the range 0.0047 ≤ q_PAH ≤ 0.047. Draine et al. (2007) found that the SINGS galaxies span the full range of q_PAH models computed, with a median value of q_PAH = 0.034. Models are further extrapolated into a (uniformly sampled) q_PAH grid, using δq_PAH = 0.001 intervals in the range 0 ≤ q_PAH ≤ 0.10, as described by AD12.

Each dust grain is assumed to be heated by radiation with an energy density per unit frequency (1)where U is a dimensionless scaling factor and is the ISRF estimated by Mathis et al. (1983) for the solar neighbourhood. A fraction (1 − γ) of the dust mass is assumed to be heated by starlight with a single intensity U = U_min, and the remaining fraction γ of the dust mass is exposed to a power-law distribution of starlight intensities between U_min and U_max, with dM/ dU ∝ U^{− α}. From now on, we call these the diffuse cloud and photodissociation regions (PDR) components respectively. AD12 found that the observed SEDs in the NGC 628 and NGC 6946 galaxies are consistent with DL models with U_max = 10⁷. Given the limited number of photometric constraints, we fix the values of U_max = 10⁷ and α = 2 to typical values found in AD12. The DL models presented in DL07 are further interpolated into a (finely sampled) U_min grid using δU_min = 0.01 intervals, as described by A15.

Therefore, in the present work the DL parameter grid has only four dimensions, Σ_Md,q_PAH,U_min, and γ. We explore the ranges 0.00 ≤ q_PAH ≤ 0.10, 0.01 ≤ U_min ≤ 30, and 0 ≤ γ ≤ 1.0. For this range of parameters, we build a DL model library that contains the model SED in a finely-spaced wavelength grid for 1 μm <λ< 1cm.

As a derived parameter, we define the ratio (2)where L_PDR is the luminosity radiated by dust in regions where U> 10² and L_dust is the total power radiated by the dust. Clearly, f_PDR depends on the fitting parameter γ in the numerator and, through the denominator, also depends on U_min. Dust heated with U> 10² emits predominantly in the λ< 100 μm range; therefore, the IRAS 60 to IRAS 100 intensity ratio can be increased to very high values by taking f_PDR → 1. Conversely, for a given U_min, the minimum IRAS 60/IRAS 100 intensity ratio corresponds to models with f_PDR = 0.

Another derived quantity, the mass-weighted mean starlight heating intensity ⟨ U ⟩, for α = 2, is given by (3)Adopting the updated carbonaceous and astrosilicate densities recommended by DA14, the DL model used here is consistent with the MW ratio of visual extinction to H column density, A_V/N_H = 5.34 × 10^-22 mag cm² (i.e. N_H/E(B − V) = 5.8 × 10²¹ cm^-2 mag^-1, Bohlin et al. 1978), for a dust to H mass ratio Σ_Md/N_Hm_H = 0.0091. From the dust surface density, we infer the model estimate of the visual extinction (4)

3.2. Fitting strategy and implementation

For each individual pixel, we find the DL parameters { Σ_Md,q_PAH,U_min,γ } that minimize (5)where S_obs(λ_k) is the observed flux density per pixel, S_DL(λ_k) is the DL emission SED convolved with the response function of the spectral band k, and σ_λk is the 1σ uncertainty in the measured intensity density at wavelength λ_k. We use a strategy similar to that of AD12 and define σ_λk as a sum in quadrature of five uncertainty sources:

• the calibration uncertainty (proportional to the observed intensity);

• the zero-level (offset) uncertainty;

• the residual dipole uncertainty;

• CIB anisotropies;

• the instrumental noise.

Values for these uncertainties (except the noise) are given in Table 1. To produce the best-fit parameter estimates, we fit the DL model to each pixel independently of the others.

We observe that for a given set of parameters { q_PAH,U_min }, the model emission is bi-linear in { Σ_Md,γ }. This allows us to easily calculate the best-fit values of { Σ_Md,γ } for a given parameter set { q_PAH,U_min }. Therefore, when looking for the best-fit model in the full four-dimensional model parameter space { Σ_Md,q_PAH,U_min,γ }, we only need to perform a search over the two-dimensional subspace spanned by { q_PAH,U_min }. The DL model emission convolved with the instrumental bandpasses, S_DL(λ_k), was pre-computed for a { q_PAH,U_min } parameter grid, allowing the multi-dimensional search for optimal parameters to be performed quickly by brute force, without relying on non-linear minimization algorithms.

In order to determine the uncertainties on the estimated parameters in each pixel, we proceed as follows: we simulate 100 observations by adding noise to the observed data; we fit each simulated SED using the same fitting technique as for the observed SED; and we study the statistics of the fitted parameters for the various realizations. The noise added in each pixel is a sum of the five contributions listed in the previous paragraph, each one assumed to be Gaussian distributed. We follow a strategy similar to that of AD12, taking a pixel-to-pixel independent contribution for the data noise and correlated contributions across the different pixels for the other four sources of uncertainty. For simplicity, we assume that none of the uncertainties are correlated across the bands. The parameter error estimate at a given pixel is the standard deviation of the parameter values obtained for the simulated SEDs. For typical pixels, the uncertainty on the estimated parameters is a few percent of their values (e.g. Fig. 2 shows the signal-to-noise (S/N) ratio of Σ_Md).

4. Dust modelling results and fitting robustness analysis

We present the results of the model fits (Sect. 4.1) and residual maps that quantify the model ability to fit the data (Sect. 4.2). In Sect. 4.1, we assess the robustness of the dust mass surface density with respect to the choice of frequency channels used in the fit.

4.1. Parameter maps

Figure 1 shows the all-sky maps of the fitted dust parameters. The left column corresponds to a Mollweide projection of the sky in Galactic coordinates, and the centre and right columns correspond to orthographic projections of the southern and northern hemispheres, centred on the corresponding Galactic poles. The top row corresponds to the dust mass surface density, Σ_Md, the main output of the model on which we focus our analysis in the next sections of the paper. Away from the Galactic plane, this map displays the structure of molecular clouds and the diffuse ISM in the solar neighbourhood. The middle row shows the map of U_min computed at the 5′ resolution of the IRAS and Planck data. At high Galactic latitude, the CIB anisotropies induce a significant scatter in U_min. Extragalactic point sources also contribute to the scatter of U_min where the Galactic dust emission is low. At low Galactic latitudes, the U_min values tend to be high (U_min> 1) in the inner Galactic disk and low (U_min< 1) in the outer galactic disk. The U_min map present structures aligned with the ecliptic plane, especially at high Galactic and low ecliptic latitudes, which are likely to be artefacts reflecting uncertainties in the subtraction of the zodiacal emission.

The f_PDR map shows artefact structures aligned with the ecliptic plane especially at high Galactic and low ecliptic latitudes. These artefacts are likely to be caused by residual zodiacal light in the IRAS 60 maps. As shown in Sect. 4.3.1, the dust mass estimates are not strongly biased in these regions. Figure 1 does not display the q_PAH maps, which are presented by A15 together with the corrected WISE data. The mass fraction in the PAH grains is relatively small, and therefore, variations in q_PAH do not have a major impact on the Σ_Md. If instead of using the WISE data to constrain q_PAH, we simply fix q_PAH = 0.04, the Σ_Md estimates will only change by a few percent.

Figure 2 shows a map of the dust emitted luminosity surface density, Σ_Ld, the mean intensity heating the dust, ⟨ U ⟩, the χ² per degree of freedom (d.o.f.) of the fit, χ²/ N.d.o.f., and a map of the S/N ratio of the dust mass surface density Σ_Md.

The χ²/ N.d.o.f. map scatter around unity in the high Galactic latitude areas, where the data uncertainties are noise-dominated. The χ²/ N.d.o.f. is slightly larger than 1 in the inner Galactic disk and several other localized areas. In the outer Galactic disk the χ²/ N.d.o.f. is smaller than 1, presumably due to overestimation of the uncertainties. Over much of the sky, the fit to the FIR SED is not as good as in Pl-MBB; the MBB fit has three fitting parameters in contrast with the DL model which has only two, Σ_Md and U_min⁹.

4.2. Dust model photometric performance: residual maps

As shown in the χ²/ N.d.o.f. map in Fig. 2, the DL model fits the observed SED satisfactorily (within 1σ) over most of the sky areas. However, the model SEDs have systematic departures from the observed SED in the inner Galactic disk, at low ecliptic latitude, and in localized regions. We note that the spectral index of the dust FIR−submm opacity is fixed in the DL model; it cannot be adjusted to match the observed SED closely. This is why MBB spectra (with one extra effective degree of freedom) fits the observed SED better in some regions. The departures of the model in the low ecliptic latitude regions could be caused by defects in the zodiacal light estimation (and removal) from the photometric maps that the model cannot accommodate. In the Magellanic Clouds (MC) the DL model fails to fit the data¹⁰. The MC exhibit surprisingly strong emission at submm and millimetre wavelengths. Planck Collaboration XVII (2011) conclude that conventional dust models cannot account for the observed 600−3000 μm emission without invoking unphysically large amounts of very cold dust. Draine & Hensley (2012) suggest that magnetic dipole emission from magnetic grain materials could account for the unusually strong submm emission from the Small MC.

Figures 3 and 4 show the model departures from the photometric constraints used in the fits. Each panel shows the difference between the model predicted intensity and the observed intensity, divided by the observed uncertainty. The systematic departures show that the physical model being used does not have sufficient parameters or flexibility to fit the data perfectly.

By increasing γ (i.e. the PDR component), the DL model can increase the IRAS 60 to IRAS 100 ratio to high values, without contributing much to the Planck intensities. Thus, in principle, the model should never underpredict the IRAS 60 emission. Figure 3 shows the model performance for fitting the IRAS bands; several high latitude areas (mostly with f_PDR = 0) have IRAS 60 overpredicted and IRAS 100 underpredicted. Both model components (the diffuse cloud and PDR components) have an IRAS 60/IRAS 100 intensity ratio slightly larger than the ratio observed in these regions. There are several areas where the IRAS 60/IRAS 100 ratio is below the value for the best-fit U_min, hence in these areas the model (with f_PDR = 0) overpredicts IRAS 60. This systematic effect is at the 1−2σ level (i.e. 10−20%).

In the inner Galactic disk the DL model tends to underpredict the 350 μm and overpredict the 850 μm emission (see Fig. 4). The observed SED is systematically steeper than the DL SED in the 350−850 μm range (i.e. between Planck 857 and Planck 353). Similar results were found in the central kiloparsec of M31 in the 250−500 μm range (DA14). The MBB fit of these regions, presented in Pl-MBB, finds larger values of the opacity spectral index β (β ≈ 2.2) than the typical value found in the low-and mid-range dust surface density areas (β ≈ 1.65). The DL SED peak can be broadened by increasing the PDR component (i.e. by raising γ or f_PDR), but it cannot be made steeper than the γ = 0 (f_PDR = 0) models, and the model therefore fails to fit the 350−850 μm SED in these regions.

Following DA14, we define (6)as the effective power-law index of the DL dust opacity between 350 μm and 850 μm, where κ_DL ∗ F is the assumed absorption cross-section per unit dust mass convolved with the respective Planck filter. For the DL model¹¹ this ratio is Υ_DL ≈ 1.8.

If the dust temperatures in the fitted DL model were left unchanged, then the predicted Planck 857/Planck 353 intensity ratio could be brought into agreement with observations if Υ_DL were changed by δΥ given by (7)where we denote I_ν(Planck...) the observed Planck intensity, and I_ν(DL ...) the corresponding intensity for the DL model.

Figure 5 shows the δΥ map, i.e. the correction to the spectral index of the submm dust opacity that would bring the DL SED into agreement with the observed SED if the dust temperature distribution is left unchanged. The observed SED is steeper than the DL model in the inner Galactic disk (δΥ_DL ≈ 0.3) and shallower in the MC (δΥ_DL ≈ −0.3). The correction to the spectral index δΥ is positive on average. The average value of δΥ tends to increase with Σ_Md, but the scatter of the individual pixels is always larger than the mean. The large dispersion in the low surface brightness areas is mainly due to CIB anisotropies. The dispersion in bright sky areas, e.g. along the Galactic plane and in molecular clouds off the plane, may be an indicator of dust evolution, i.e. variations in the FIR emission properties of the dust grains in the diffuse ISM.

Modifying the spectral index of the dust opacity in the model would change U_min and thereby the dust mass surface density. The δΥ map should be regarded as a guide on how to modify the dust opacity in future dust models, rather than as the exact correction to be applied to the opacity law per se.

4.3. Robustness of the mass estimate

4.3.1. Importance of IRAS 60

To study the potential bias introduced by IRAS 60, due to residuals of zodiacal light estimation (whose relative contribution is the largest in the IRAS 60 band) or the inability of the DL model to reproduce the correct SED in this range, one can perform modelling without the IRAS 60 constraint. In this case we set γ = 0, i.e. we allow only the diffuse cloud component (f_PDR = 0), and so we have a two-parameter model.

Figure 6 shows the ratio of the dust mass estimated without using the IRAS 60 constraint and with γ = 0 to that estimated using IRAS 60 and allowing γ to be fitted (i.e. our original modelling). The left panel shows all the sky pixels and the right panel only the pixels with f_PDR> 0. In the mid-and-high-range surface mass density areas (Σ_Md> 10⁵M_⊙ kpc^-2), where the photometry has good S/N, both models agree well, with a rms scatter below 5%. The inclusion or exclusion of the IRAS 60 constraint does not significantly affect our dust mass estimates in these regions. In the low surface density areas, inclusion of the IRAS 60 does not change the Σ_Md estimate in the f_PDR> 0 areas, but it leads to an increase of the Σ_Md estimate in the f_PDR = 0 pixels. In the f_PDR = 0 areas, the model can overpredict IRAS 60 in some pixels, and therefore, when this constraint is removed, the dust can be fitted with a larger U_min value reducing the Σ_Md needed to reproduce the remaining photometric constraints. In the f_PDR> 0 areas, the PDR component has a small contribution to the longer wavelengths constraints, and therefore removing the IRAS 60 constraint and PDR component has little effect in the Σ_Md estimates.

4.3.2. Dependence of the mass estimate on the photometric constraints

The Planck and IRAS data do not provide photometric constraints in the 120 μm <λ< 300 μm range. This is a potentially problematic situation, since the dust SED typically peaks in this wavelength range. We can add the DIRBE 140 and DIRBE 240 constraints in a low resolution (FWHM> 42′) modelling to test this possibility.

We compare two analyses performed using a 1° FWHM Gaussian PSF. The first uses the same photometric constraints as the high resolution modelling (WISE, IRAS, and Planck), and the second additionally uses the DIRBE 140 and DIRBE 240 constraints. The results are shown in Fig. 7. Both model fits agree very well, with differences between the dust mass estimates of only a few percent. Therefore, our dust mass estimates are not substantially affected by the lack of photometric constraints near the SED peak. This is in agreement with similar tests carried out in Pl-MBB.

4.3.3. Validation on M31

In Appendix A we compare our dust mass estimates in the Andromeda galaxy (M31) with estimates based on an independent data set and processing pipelines. Both analyses use the DL model. This comparison allows us to analyse the impact of the photometric data used in the dust modelling. We conclude that the model results are not sensitive to the specific data sets used to constrain the FIR dust emission, validating the present modelling pipeline and methodology.

5. Comparison between the DL and MBB optical exctinction estimates

We now compare the DL optical extinction estimates with those from the MBB dust modelling presented in Pl-MBB, denotedA_V,DL andA_V,MBB respectively. We perform a DL dust modeling of the same Planck data as in Pl-MBB (i.e. the nominal mission maps)¹², the same IRAS 100 data, but our DL modelling also includes IRAS 60 and WISE 12 constraints. The DL model has two extra parameters (γ and q_PAH) that can adjust the IRAS 60 and WISE 12 intensity fairly independently of the remaining bands. Therefore, the relevant data that both models are using in determining the FIR emission are essentially the same. The MBB extinction map has been calibrated with external (optical) observational data, and so this comparison allows us to test the DL modelling against those independent data.

Pl-MBB estimated the optical extinction¹³A_V,QSO for a sample of QSOs from the SDSS survey. A single normalization factor Π was chosen to convert their optical depth τ₃₅₃ map (the parameter of the MBB that scales linearly with the total dust emission, similar to the DL Σ_Md) into an optical extinction map: A_V,MBB ≡ Π τ₃₅₃.

DL is a physical dust model and therefore fitting the observed FIR emission directly provides an optical extinction estimate, without the need for an extra calibration. However, if the DL dust model employs incorrect physical assumptions (e.g. the value of the FIR opacity), it may systematically over or under estimate the optical extinction corresponding to observed FIR emission.

Figure 8 shows the ratio of the DL and MBB A_V estimates. The top row shows the ratio map. The bottom row shows its scatter and histogram. Over most of the sky (0.1 mag <A_V,DL< 20 mag), the A_V,DL values are larger than the A_V,MBB by a factor of 2.40 ± 0.40. This discrepancy is roughly independent of A_V,DL. The situation changes in the very dense areas (inner Galactic disk). In these areas (A_V,DL ≈ 100 mag), the A_V,DL are larger than the A_V,MBB estimates by 1.95 ± 0.10.

In the diffuse areas (A_V,DL ≲ 1), where the A_V,MBB has been calibrated using the QSOs, A_V,DL overestimates A_V,MBB, and therefore A_V,DL should overestimate A_V,QSO by a similar factor. This mismatch arises from two factors. (1) The DL dust physical parameters were chosen so that the model reproduces the SED proposed by Finkbeiner et al. (1999), based on FIRAS observations. It was tailored to fit the high latitude I_ν/N(H i) with U_min ≈ 1. The high latitude SED from Planck observations differs from that derived from FIRAS observations. The difference depends on the frequency and can be as high as 20% (Planck Collaboration VIII 2014). The best fit to the mean Planck + IRAS SED on the QSO lines of sight is obtained for U_min ≈ 0.66. The dust total emission (luminosity) computed for the Planck and for the Finkbeiner et al. (1999) SED are similar. The dust total emission per unit of optical reddening (or mass) scales linearly with U_min. Therefore, we need 1/0.66 ≈ 1.5 more dust mass to reproduce the observed luminosity. This is in agreement with the results of Planck Collaboration Int. XVII (2014) who have used the dust −H i correlation at high Galactic latitudes to measure the dust SED per unit of H i column density. They find that their SED is well fit by the DL model for U_min = 0.7 after scaling by a factor 1.45. (2) The optical extinction per gas column density used to construct the DL model is that of Bohlin et al. (1978). Recent observations show that this ratio needs to be decreased by a factor of approximately 1/1.4 (Liszt 2014a,b). Combining the two factors, we expect the A_V,DL to overestimate the A_V,QSO by about a factor of 2.

6. Renormalization of the model extinction map

We proceed in our analysis of the model results characterizing how the ratio between the optical extinction from the DL model and that measured from the optical photometry of QSOs depends on U_min. We introduce the sample of QSOs we use in Sect. 6.1, and compare the DL and QSO A_V estimates in Sect. 6.2. Based on this analysis, we propose a renormalization of the optical extinction derived from the DL model (Sect. 6.3). The renormalized extinction map is compared to that derived by Schlafly et al. (2014, hereafterSGF) from stellar reddening using the Pan-STARRS1 (Kaiser et al. 2010) data.

6.1. The QSO sample

SDSS provides photometric observations for a sample of 272 366 QSOs, which allow us to measure the optical extinction for comparison with that from the DL model. A subsample of 105 783 (an earlier data release) was used in Pl-MBB to normalize the opacity maps derived from the MBB fits in order to produce an extinction map.

The use of QSOs as calibrators has several advantages over other cross-calibrations:

• QSOs are extragalactic, and at high redshift, so all the detected dust in a given pixel is between the QSO and us, a major advantage with respect to maps generated from stellar reddening studies;

• the QSO sample is large and well distributed across diffuse (A_V ≲ 1) regions at high Galactic latitude, providing good statistics;

• SDSS photometry is very accurate and well understood.

In Appendix B we describe the SDSS QSO catalogue in detail, and how for each QSO we measure the extinction A_V,QSO from the optical SDSS observations.

6.2. A_V,QSO – A_V,DL comparison

In this section, we present a comparison of the DL and QSO extinction, as a function of the fitted parameter U_min. The DL and QSO estimates of A_V are compared in the following way. We sort the QSO lines of sight with respect of the U_min value and divide them in 20 groups having (approximately) equal number of QSOs each. For each group of QSOs, we measure the slope ϵ(U_min) fitting the A_V,QSO versus A_V,DL data with a line going through the origin. In Fig. 9, ϵ(U_min) is plotted versus the mean value of U_min for each group. We observe that ϵ(U_min), a weighted mean of A_V,QSO/A_V,DL in each U_min bin, depends on U_min. The slope obtained fitting the A_V,QSO versus A_V,DL for the whole sample of QSOs is ⟨ ϵ ⟩ ≈ 0.495. Therefore, on average the DL model overpredicts the observed A_V,QSO by a factor of 1/0.495 = 2.0, with the discrepancy being larger for sightlines with smaller U_min values. There is a 20% difference between the 2.4 factor that arises from the comparison between A_V,DL and A_V,QSO indirectly via the MBB A_V fit, and the mean factor of 2.0 found here. This is due to the use of a different QSO sample (Pl-MBB used a smaller QSO sample), which accounts for 10% of the difference, and to differences in the way the QSO A_V is computed from the SDSS photometry, responsible of the remaining 10%.

For a given FIR SED, the DL model predicts the optical reddening unambiguously, with no freedom for any extra calibration. However, if one had the option to adjust the DL extinction estimates by multiplying them by a single factor (i.e. ignoring the dependence of ϵ on U_min), one would reduce the optical extinction estimates by a factor of 2.0.

6.3. Renormalization of the A_V,DL map in the diffuse ISM

We use the results of the A_V,QSO –A_V,DL comparison (Sect. 6.2 and Fig. 9) to renormalize the A_V,DL map in the diffuse ISM. The ratio between the two extinction values is well approximated as a linear function of U_min: (8)Thus, we define a renormalized DL optical reddening as¹⁴: (9)Empirically, A_V,RQ is our best estimator of the QSO extinction A_V,QSO. Pl-MBB proposed the dust radiance (the total luminosity emitted by the dust) as a tracer of dust column density in the diffuse ISM. This would be expected if the radiation field heating the dust were uniform, and the variations of the dust temperature were only driven by variation of the dust FIR-submm opacity in the diffuse ISM. The dust radiance is proportional to U_min × A_V,DL. Our best fit of the renormalization factor as a function of U_min is an intermediate solution between the radiance and the model column density (A_V,DL). The scaling factor in Eq. (8) increases with U_min but with a slope smaller than 1. Figure 9 shows that our renormalization is a better fit of the data than the radiance.

6.4. Validation of the renormalized extinction with stellar observations

We compare the renormalized extinction A_V,RQ to the extinction estimated by SGF from optical stellar observations, denoted A_V,Sch. SGF presented a map of the dust reddening to 4.5 kpc derived from Pan-STARRS1 stellar photometry. Their map covers almost the entire sky north of declination −30° at a resolution of 7″−14″. In the present analysis, we discard the sky areas with | b | < 5° to avoid the Galactic disk where a fraction of the dust is farther than 4.5 kpc, and therefore, traced by dust in emission, but not present in the SGF absorption map.

Figure 10 shows the comparison between the renormalized DL A_V estimates (A_V,RQ) with the stellar observations based A_V,Sch estimates in the | b | > 5° sky. The agreement between both A_Vestimates is remarkable. The difference between the estimates (A_V,Sch−A_V,RQ) has a mean of 0.02 mag and scatter (variance) of 0.12 mag. This comparison validates our renormalized DL A_V estimates as a good tracer of the dust optical extinction in the | b | > 5° sky. Our optical extinction estimate (A_V,RQ) was tailored to match the QSOs A_V estimates. QSOs were observed in diffuse diffuse Galactic lines of sight (most of the QSOs have A_V ≈ 0.1). The agreement of the A_V,Sch and A_V,RQ estimates extends the validity of the DL renormalized estimate (A_V,RQ) to greater column densities. The scatter of A_V,Sch–A_V,RQ provides an estimate of the uncertainties in the different A_V estimates.

7. FIR SEDs per unit of optical extinction

The DL model parameter U_min and the QSO data are used to compress the FIR IRAS and sub-mm Planck observations of the diffuse ISM to a set of 20 SEDs normalized per A_V, which we present and discuss.

The parameter U_min is mainly determined by the wavelength where the SED peaks; as a corollary, SEDs for different values of U_min differ significantly. The A_V values obtained from the QSO analysis,A_V,QSO, allow us to normalize the observed SEDs (per unit of optical extinction) and generate a one-parameter family of I_ν/A_V. This family is indexed by the U_min parameter; the QSO lines of sight are grouped according to the fitted Galactic U_min value. We divide the sample of good QSOs in 20 bins, containing 11 212 QSOs each¹⁵.

To obtain the I_ν/A_V values we proceed as follows. For each band and U_min, we would like to perform a linear regression of the I_ν values as a function ofA_V,QSO. The large scatter and non-Gaussian distribution ofA_V,QSO and the scatter on I_ν make it challenging to determine such a slope robustly. Therefore, we smooth the maps to a Gaussian PSF with 30′ FWHM to reduce the scatter on I_ν, redo the dust modelling (to obtain a coherent U_min estimate), and perform the regression on the smoothed (less noisy) maps. The non-Gaussian distribution ofA_V,QSO do not introduce any bias in the slope found¹⁶.

Figure F11 presents the set of SEDs. The left panel shows the SEDs for the different U_min values. The right panel shows each SED divided by the mean to highlight the differences between the individual SEDs.

The complex statistics of theA_V,QSO estimates, which depend on variations in QSO intrinsic spectra, makes it hard to obtain a reliable estimate of the uncertainties in the (normalizing) A_V. However, the statistical uncertainties are not dominant because they average out thanks to the large size of the QSOs sample. The accuracy of our determination of the FIR intensities per unit of optical extinction is mainly limited by systematic uncertainties on the normalization by A_V. Pl-MBB estimated the A_V over a subsample of the QSOs, and their estimates differ from ours by ≈14%. When we compare ourA_V,QSO estimates from those of SGF based on stellar photometry (Sect. 6.4), we find an agreement within 5−10% over the QSO lines of sight, Therefore, it is reasonable to assume ourA_V,QSO estimates are uncertain at a 10% level. The instrumental calibration uncertainties (≲6% at Planck frequencies) translate directly to the FIR intensities per unit of optical extinction. Therefore, the normalization of each SED may be uncertain up to about ≈15%.

In Fig. 12, the specific intensities per unit of optical extinction are compared with the DL model. The four panels correspond to the spectral bands: IRAS 100; Planck 857; Planck 545; and Planck 353. In each panel, the black curve corresponds to the DL intensity normalized by A_V, and the red curves to the DL intensity normalized by A_V,RQ, i.e. the model intensity scaled by the renormalization factor in Eq. (8). The DL model under-predicts the FIR intensities per unit of optical extinction A_V by significant amounts, especially for sightlines with low fitted values of U_min, but the renormalization of the DL A_V values brings into agreement the observed and model band intensities per unit of extinction. This result shows that the DL model has approximately the correct SED shape to fit the measured SEDs for the diffuse ISM, and that a U_min-dependent renormalization brings the DL model into agreement with the IRAS and Planck data. It is also a satisfactory check of the consistency of our data analysis and model fitting. Indeed, the SED values in Fig. 12 are derived from a linear fit between the data andA_V,QSO, the renormalization factor from a linear fit betweenA_V,QSO andA_V,DL (see Fig. 9), andA_V,DL from the DL model fit of the data (see Sect. 3.2).

8. Optical extinction and FIR emission in molecular clouds

We extend our assessment of the DL extinction maps to molecular clouds using extinction maps from 2MASS stellar colours that are presented in Sect. 8.1. In Sect. 8.2, we compare them with the original and renormalized estimates derived from the DL model. We discuss a model renormalization for molecular clouds in Sect. 8.3.

8.1. Extinction maps of molecular clouds

Schneider et al. (2011) presented optical extinction maps, denoted A_V,2M, of several clouds computed using stellar observations from the 2MASS catalogue in the J, H, and K bands. We use the maps of the Cepheus, Chamaeleon, Ophiuchus, Orion, and Taurus cloud complexes. The Schneider et al. (2011) A_V maps were computed using a 2′ Gaussian PSF, and we smooth them to a 5′ Gaussian PSF to perform our analysis. We corrected the 2MASS maps for a zero level offset that is fitted with an inclined plane with an algorithm similar to the one used to estimate the background in the analysis of the KINGFISH sample of galaxies in AD12. The algorithm iteratively and simultaneously matches the zero level across the A_V,RQ and A_V,2M maps and defines the areas that are considered background. The uncertainty on the zero level of the A_V maps is smaller than 0.1 magnitude. It is significant only for the map areas with the lowest A_V values.

Figure F13 shows the 2MASS A_V,2M map, the DL U_min map, theA_V,DL map (divided by 3.07, see Sect. 8.2 ) and the renormalized A_V,RQ map for the Chamaeleon region. The inner (high A_V) areas correspond to the lowest U_min values, as expected since the stellar radiation field heating dust grains is attenuated when penetrating into molecular clouds. The cloud complexes show a similar A_V −U_min trend.

8.2. Comparison of 2MASS and DL extinction maps in molecular clouds

For each cloud, we find an approximate linear relation between the A_V,2M and theA_V,DL maps as illustrated for the Chamaeleon cloud in the left panel of Fig. 14. After multiplicative adjustment, theA_V,DL and A_V,2M estimates agree reasonably well. However, as in the diffuse ISM, the (FIR based)A_V,DL estimates are significantly larger than the (optical) A_V,2M estimates. For the selected clouds, the DL model overestimates the 2MASS stellar A_V by factors of 2−3. Table 2 provides the multiplicative factors needed to make the 2MASS A_V maps agree with the DL A_V maps.

We compare the renormalized A_V,RQ versus A_V,2M values in the right panel of Fig. 14 for the Chamaeleon cloud and in Fig. 15 for all of the clouds. We find that the model renormalization, computed to bring into agreement the DL and QSO A_V estimates in the diffuse ISM, also accounts quite well (within 10%) for the discrepancies between 2MASS and DL A_V estimates in molecular clouds in the 0 <A_V< 3 range, and even does passably well (within 30%) up to A_V ≈ 8.

8.3. Renormalization ofA_V,DL in molecular clouds.

In the diffuse ISM analysis, we concluded that U_min is tracing variations in the radiation field and the dust opacity. Both phenomena are also present in molecular clouds, but their relative contribution in determining the SED peak (and therefore U_min) need not be the same as in the diffuse ISM. Therefore, a renormalization of the DL A_V based on 2MASS data may be different from that determined using the QSOs.

In Fig. 16 we compare the DL and 2MASS A_V estimates for pixels from the five cloud complexes with 1 <A_V,2M< 5. Similarly to Fig. 9, but using the 2MASS A_V estimates instead of those from QSOs, we plot the ratioA_V,DL/A_V,2M as a function of the fitted U_min. For each U_min value, the solid curve corresponds to the best fit slope of theA_V,DL versus A_V,2M values (i.e. it is an estimate of a weighted mean of theA_V,DL/A_V,2Mratio). The straight solid line corresponds to a fit to the solid curve in the 0.2 <U_min< 1.0 range. In this fit, each U_min is given a weight proportional to the number of pixels that have this value in the clouds (i.e. most of the weight is for the pixels within the 0.2 <U_min< 0.8 range). The dashed line corresponds to the renormalization proposed in Sect. 6.3 (Eq. (9)) for the diffuse ISM.

The straight line in Fig. 16 corresponds to a renormalization tailored to bring into agreement theA_V,DL and A_V,2M estimates, i.e. a 2MASS renormalization for molecular clouds, denoted A_V,RC. The 2MASS renormalization is given by (10)Empirically, A_V,RC is our best estimator of the 2MASS extinction A_V,2M. It is satisfactory for the renormalization method to find that the 2MASS normalization factor for molecular clouds is quite close to the one for the diffuse ISM.

9. Discussion

In Sects. 9.1 and 9.2, we relate the difference between the DL model and A_V estimates from QSO observations to dust emission properties and their evolution within the diffuse ISM. The renormalization method we propose to correct empirically for this discrepancy is discussed in Sect. 9.3.

9.1. DL FIR emission and optical extinction disagreement

In the diffuse ISM, the DL model provides good fits to the SED of Galactic dust from WISE, IRAS, and Planck data, as it has been the case in the past for external galaxies observed with Spitzer and Herschel. However, the fit is not fully satisfactory because the A_V values from the model do not agree with those derived from optical colours of QSOs for the diffuse ISM and stars for molecular clouds. The optical extinction discrepancy can be decomposed in two levels: (1) the mean factor of 1.9 between the DL and QSOs A_V values; and (2) the dependence of the ratio between the DL and QSO A_V values on U_min.

The result of the SED fit depends on the spectral shape of the dust opacity. A spectral difference makes the DL model fit with a lower U_min value than the true radiation field intensity, which turns into an increase of the A_V estimates¹⁷. The mean factor between the DL and QSOs A_V values could also be indicating that the DL dust material has a FIR–submm opacity per unit of optical extinction that is too low.

The dependency of the ratio between the DL and QSO A_V values on U_min shows that variations in the value of the FIR−submm dust opacity per unit of optical extinction, or its spectral shape, are needed across the sky; we take this to be evidence of dust evolution. This discrepancy will be present for all dust models based on fixed dust optical properties, possibly with a different magnitude depending of the details of the specific model.

Fanciullo et al. (2015) have used three dust models, the DL, Compiègne et al. (2011) and Jones et al. (2013) models, to fit the Planck SEDs in Table 3 and compare results. We note that the SEDs listed in the Table and used by Fanciullo et al. (2015) are the DL model results obtained fitting the same Planck and IRAS data as in Pl-MBB (see Sect. 2). This study, a follow-up of our data analysis, confirms our interpretation of our DL modelling in two ways. First, the mean ratio between the model and QSO A_V values is different for each model. The best match between model and data is obtained for the model of Jones et al. (2013), which uses a FIR–submm opacity per unit of optical extinction¹⁸ larger than that of the DL model. Second, the three models fail in the same way to reproduce the variations in the emission per unit extinction with a variable ISRF intensity.

Fanciullo et al. (2015) also present fits of the Planck SEDs in Table 3 with a MBB spectrum (see their Appendix A). The dust FIR-submm opacity per unit extinction and the temperature obtained from these empirical fits display the same anti-correlation as that observed for the dust models between A_V and the ISRF intensity. Thus, the dust models and MBB fits provide the same evidence for changes in the FIR-submm dust opacity in the diffuse ISM. Fanciullo et al. (2015) estimate that the amplitude of the variations must be ~20% to account for the typical (± 1σ around the mean) variations of the dust SEDs, and 40−50% for the full range.

9.2. Dust evolution in the diffuse ISM

We relate the Planck evidence for variations of the FIR–submm dust opacity to earlier studies of dust evolution in the ISM.

The extinction curve is known to vary through the ISM (Cardelli et al. 1989; Fitzpatrick & Massa 2007), especially in the UV but also in the near-IR (Fitzpatrick & Massa 2009). Similarly, the mid-IR dust emission presents variations of colour ratios between bands (Miville-Deschênes et al. 2002; Flagey et al. 2009). It is one of the successes of models like DL07 to be able accommodate this kind of evolution by adapting the size distribution of each component (Weingartner & Draine 2001), as well as the radiation field intensity for emission.

With the observations of FIR dust emission by the ISOPHOT camera onboard ISO and the PRONAOS balloon, it was however demonstrated that adapting the size distribution could not explain the low grain temperatures observed in dense filaments (Del Burgo et al. 2003; Stepnik et al. 2003). These papers argue that the intrinsic dust opacity in the FIR increases by a factor of a few, and that this can be achieved by modifying the grain structure and composition itself. Indeed, dust is known to be more emissive as composite aggregates than as compact grains (Stognienko et al. 1995; Köhler et al. 2011). The number of observations has since increased and confirmed this tendency (Del Burgo & Laureijs 2005; Ridderstad et al. 2006; Planck Collaboration XXV 2011), but primarily in the dense ISM (A_V larger than a few magnitudes). More recently with Spitzer and Planck observations such variations in the dust FIR opacity have been identified in the diffuse ISM, with amplitudes comparable to those observed before in the dense medium (Bot et al. 2009; Martin et al. 2012; Planck Collaboration XI 2014; Planck Collaboration Int. XVII 2014).

It is reasonable to accommodate the increase of dust opacity in the dense ISM by invoking coagulation (Köhler et al. 2011; Ysard et al. 2012, 2013). However, dust models face now a new challenge. Since the timescale for dust coagulation increases for decreasing density, this solution does not hold in the diffuse, low-density, ISM. Variations in grain mantles composition and thickness through the accretion of carbon and cumulative irradiation have been proposed by Jones et al. (2013) and Ysard et al. (2015) as a main dust evolutionary process within the diffuse ISM. Differences in the past history of dust grains, i.e. in their evolutionary cycle between the diffuse ISM and molecular clouds, could also contribute to the observed scatter in the shape, size and composition of dust grains in the diffuse ISM (Martin et al. 2012). Modelling is needed to test whether these ideas can match the signature of dust evolution reported in this paper.

9.3. Optical exctinction A_V estimates of dust models

To obtain an accurate A_V map of the sky, we proposed two renormalizations of the A_V estimates derived from the DL model (A_V,RQ in Eq. (9) and A_V,RC in Eq. (10)) that compensates for the discrepancy between the observed FIR emission and the optical extinction for the diffuse ISM and molecular clouds. Essentially, the renormalization rescales one of the model outputs (the dust optical extinction A_V) by a function of U_min, to match data. Planck Collaboration Int. XXVIII (2015) presented an independent comparison of the renormalized A_V,RQ estimates with γ-ray observations in the Chamaeleon cloud. They concluded that the renormalized A_V,RQ estimates are in closer agreement with γ-ray A_V estimates than the (non-renormalized)A_V,DL estimates. We now discuss the model renormalization in a more general context.

The renormalized DL estimates (A_V,RQ and A_V,RC) provide a good A_V determination in the areas where they were calibrated, but they do not provide any insight into the physical dust properties per se; the renormalized dust model becomes simply a family of SEDs used to fit the data, from which we construct and calibrate an observable quantity (A_V,RQ and A_V,RC). Unfortunately, the fitted parameters of the renormalized model (U_min) lack a physical interpretation: U_min is not solely tracing the heating intensity of the radiation field, as assumed in DL.

The A_V estimate of the DL dust model is a function of its fitted parameters, i.e. A_V = f(Σ_Md,q_PAH,γ,U_min). In general, if we fit a dust model with several parameters, A_V will be a function of the most relevant parameters¹⁹. The DL model assumes A_V = f(Σ_Md) = k × Σ_Md, with . Our proposed renormalizations are a first step towards a functional renormalization by extending A_V = k × Σ_Md into A_V = g(U_min) × k × Σ_Md, where we take g(U_min) to be a linear function of U_min. Due to the larger scatter in the QSO A_V estimates, only a simple linear function g(U_min) can be robustly estimated in the diffuse ISM. In molecular clouds, where the data are less noisy, one could find a smooth function g′(U_min), which better matches theA_V,DL/A_V,2M fit for each U_min (i.e. in Fig. 16, the solid curve flattens for U_min> 0.8, departing from its linear fit).

Unfortunately any renormalization procedure, while leading to a more accurate A_V estimate, does not provide any further insight into the dust physical properties. Real physical knowledge will arise from a new generation of dust models that should be able to predict the correct optical extinction A_V from first principles. The next generation of dust models should be able to fit the empirical SEDs presented in Sect. 7 directly. While such a new generation of dust models is not yet available, we can for now correct for the systematic departures via Eqs. (9) and (10) for the diffuse ISM and molecular clouds, respectively.

10. Conclusions

We present a full-sky dust modelling of the new Planck data, combined with ancillary IRAS and WISE data, using the DL dust model. We test the model by comparing these maps with independent estimates of the dust optical extinction A_V using SDSS QSO photometry and 2MASS stellar data. Our analysis provides new insight on interstellar dust and a new A_V map over the full sky.

The DL model fits the observed Planck, IRAS, and WISE SEDs well over most of the sky. The modelling is robust against changes in the angular resolution, as well as adding DIRBE 140 and DIRBE 240 photometric constraints. The high resolution parameter maps that we generated trace the Galactic dusty structures well, using a state-of-the-art dust model.

In the diffuse ISM, the DL A_V estimates are larger than estimates from QSO optical photometry by approximately a factor of 2, and this discrepancy depends systematically on U_min. In molecular clouds, the DL A_V estimates are larger than estimates based on 2MASS stellar colours by a factor of about 3. Again, the discrepancy depends in a similar way on U_min.

We conclude that the current parameter U_min, associated with the peak wavelength of the SED, does not trace only variations in the intensity of the radiation field heating the dust; U_min also traces dust evolution: i.e. variations in the optical and FIR properties of the dust grains in the diffuse ISM. DL is a physical dust model. Physical dust models have the advantage that, if successful, they give some support to the physical assumptions made about the interstellar dust and ISM properties that they are based on. Unfortunately, the deficiency found in this study indicates that some of the physical assumptions of the model need to be revised.

We provide a one-parameter family of SEDs per unit of dust optical extinction in the diffuse ISM. These SEDs, which relate the dust emission and absorption properties, are independent of the dust/gas ratio or problems inferring total H column density from observations. The next generation of dust models will need to reproduce these new SED estimates.

We propose an empirical renormalization of the DL A_V map as a function of the DL U_min parameter. The proposed renormalization (A_V,RQ), derived to match the QSO A_V estimates for the diffuse ISM, also brings into agreement the DL A_V estimates with those derived from stellar colours for the | b | > 5 deg sky using the Pan-STARRS1 data (SGF), and towards nearby molecular clouds in the 0 <A_V< 5 range using the 2MASS survey. We propose a second renormalized DL A_V estimate (A_V,RC) tailored to trace the A_V estimates in molecular clouds more precisely.

The renormalized map A_V,RQ based on our QSOs analysis is our most accurate estimate of the optical extinction in the diffuse ISM. Comparison of the A_V,RQ map against other tracers of interstellar extinction that probe different environments, would further test its accuracy and check for any potential systematics. A comparison with Fermi data towards the Chamaeleon molecular cloud shows that the A_V,RQ map more closely matches the γ-ray diffuse emission than the 353 GHz opacity and radiance maps from Pl-MBB, but not as well as the fit obtained combining H i, CO, and dark neutral medium maps, which indicates significant differences in the FIR-submm dust emission properties between these gas components, not taken into account in our renormalization (Planck Collaboration Int. XXVIII 2015).

---

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.

References

Appendix A: Comparison with Spitzer + Herschel modelling of the Andromeda galaxy

The Andromeda galaxy is the nearest large spiral galaxy. It provides a useful benchmark to validate the current dust modelling. Its isophotal radius is R₂₅ = 95′ (de Vaucouleurs et al. 1991), corresponding to R₂₅ = 20.6 kpc at the assumed distance d = 744 kpc (Vilardell et al. 2010).

Several authors have modelled the dust properties of M31. Planck Collaboration Int. XXV (2015) presented an independent study to M31 using Planck maps and MBB dust model. In particular DA14 presented a DL based modelling of M31 using the IRAC (Fazio et al. 2004) and MIPS (Rieke et al. 2004) instruments on Spitzer, and the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments on Herschel. This data set has 13 photometric constraints (IRAC 3.6 μm, 4.5 μm, 5.8 μm, and 8.0 μm, MIPS 24 μm, 70 μm, and 160 μm, PACS 70 μm, 100 μm, and 160 μm, and SPIRE 250 μm, 350 μm, and 500 μm) from a different set of instruments than those used in our analysis. The high resolution modelling traces the structures of M31 in great detail, providing maps of U_min and dust surface density, and enables a comparison to be made with gas and metallicity observations. The modelling techniques are described and validated on NGC 628 and NGC 6946 in AD12, and later expanded to the full KINGFISH galaxy sample in AD15.

We compare the dust mass surface density maps²⁰ of the modelling presented by DA14 (from now on called “Herschel”) degraded to a 5′ Gaussian PSF, with the current modelling, called “Planck”. In the Herschel modelling, a tilted plane is fitted to the background areas, and subtracted from the original images to remove the Milky Way cirrus emission. Therefore, we need to add the cirrus emission back to the Herschel mass estimates before comparing to the Planck modelling. The zero level of the Herschel modelling was restored with an algorithm similar to that used to estimate the background planes in the KINGFISH dust modelling (see AD12). This algorithm iteratively fits an inclined plane to the difference in mass surface densities over the background points.

M31 does not have considerable quantities of cold dust, which would be detected in the Planck modelling but not in the Herschel modelling. Therefore, we expect both modellings to agree well.

Figure A.1 presents the comparison of the two dust models. The Herschel and Planck approaches agree very well: the resolved mass differences between the two analyses is small, only 10 % across most of the galaxy. The remaining parameter estimates also agree well. In conclusion, the model results appear not to be sensitive to the specific data sets used to constrain the FIR dust emission. This comparison validates the present modelling pipeline and methodology.

Appendix B: QSO A_v estimation

The intrinsic colours of an unobscured QSO depend strongly on its redshift²¹ (ζ). We first estimate the (redshift dependent) unobscured QSO colour for each band pair. By comparing each QSO colours with the expected unobscured colours, we can estimate its reddening. Assuming a typical dust extinction curve, we can combine the reddening estimates of the band pairs into a single extinction estimate for each QSO. This analysis relies on the fact that the mean colour excess of a group of QSO scales linearly with the DL A_V estimates (see Fig. B.2).

Appendix B.1: SDSS QSO catalogue

The SDSS is a photometric and spectroscopic survey, using a dedicated 2.5-m telescope at Apache Point Observatory in New Mexico. It has produced high quality observations of approximately 10⁴ deg² of the northern sky in five optical and near-IR bands: u, g, r, i, and z, centred at 354.3 nm, 477.0 nm, 623.1 nm, 762.5 nm, and 913.4 nm respectively (York et al. 2000). The SDSS seventh data release (DR7, Abazajian et al. 2009) contains a sample of 105 783 spectroscopically confirmed QSOs, and the SDSS tenth data release (DR10, Pâris et al. 2014) contains an additional sample of 166 583 QSOs.

In order to avoid absorption from the intergalactic medium, each SDSS band is only usable up to the redshift at which the Lyα line (121.57 nm vacuum wavelength) enters (from the blue side) into the filter. Therefore, we can use the u-band, for QSOs with ζ< 1.64, g-band for ζ< 2.31, r-band for ζ< 3.55, i-band for ζ< 4.62, and z-band for ζ< 5.69. We also limit the study to 0.35 <ζ< 3.35, to have enough QSOs per unit of redshift to estimate reliably the redshift-dependent unobscured QSO intrinsic colour (see Section B.2). We also remove the few QSOs that lie in a line of sight where the Galactic dust is very hot (U_min> 1.05), very cold (U_min< 0.4), very luminous (L_dust> 10⁸L_⊙ kpc^-2), or where there is strong extinction (A_V> 1). This leaves 261 841 useful QSOs.

Appendix B.2: Unobscured QSO intrinsic colours and extinction estimation

A typical QSO spectrum has several emission and absorption lines superimposed on a power-law-like continuum. Depending on the QSO redshift, the lines fall in different filters. Therefore, for each optical band pair (X,Y), the unobscured QSO intrinsic colour C_X,Y(ζ) depends on the QSO redshift. Given two photometric bands X and Y, in order to estimate the unobscured QSO intrinsic colour C_X,Y(ζ), we proceed as follows.

We will see that the intrinsic dust properties appear to depend on the parameter U_min. Therefore, to avoid introducing a potential bias when computing C_X,Y(ζ), we group the lines of sight according to U_min, and analyse each group independently. The functions C_X,Y(ζ) should, in principle, not depend on U_min, and therefore, all the estimates C_X,Y(ζ,U_min) should be similar for the different U_min sets. Working independently on each U_min, for each redshift ζ we choose all the QSOs in the interval [ ζ − 0.05,ζ + 0.05 ], or the 2000 closest QSOs if there are more than 2000 QSOs in the interval, and fit the QSOs colour (X − Y) as a function of the dust column density: (B.1)whereA_V,DL is the DL estimated dust extinction in each QSO line of sight. The function C_X,Y(ζ,U_min) is the best estimate of the colour difference (X − Y) of an unobscured QSO (A_V,DL = 0) at redshift ζ, estimated from the lines of sight of dust fitted with U_min. The function η_X,Y(ζ,U_min) should be essentially independent of ζ²². Variations in the function η_X,Y(ζ,U_min) with respect to U_min give us information about the dust properties.

Once we compute C_X,Y(ζ,U_min) for the different values of U_min, we average them for each redshift ζ to obtain C_X,Y(ζ). For each U_min and ζ, the weight given to each C_X,Y(ζ,U_min) value is proportional to the number of QSO in the [ ζ − 0.05,ζ + 0.05 ] interval. Figure B.1 shows the results of this unobscured QSO intrinsic colour estimation algorithm for the bands i and z. The functions C_i,z(ζ,U_min) are shown for the different values of U_min, using redder lines for larger U_min, and greener for smaller U_min. Their weighted mean C_i,z(ζ) is shown in black.

For each QSO, we define its reddening E_X,Y as: (B.2)The E_X,Y values should not depend on the redshift, and therefore we can group all the QSOs of a given U_min into a sub sample with the same intrinsic colour. We note that no additional hypotheses on the QSO spectral shape or dust extinction curve need to be made to compute the QSO intrinsic colours. Working with all the QSOs with a given U_min, we fit (B.3)and identify the outlier QSOs that depart by more than 3σ from the expected linear relationship. Figure B.2 shows the typical QSO E_g,r versusA_V,DL fit for U_min = 0.6. In this case, η_g,r(U_min = 0.6) = 0.19. Although the QSO E_X,Y versusA_V,DL relationship has large scatter due to variations in the QSOs spectra (continuum and lines) and intrinsic obscuration in the QSOs, as long as there is no selection bias with respect toA_V,DL our study should be robust. The fact that the mean QSO E_X,Y for eachA_V,DL (curve) and the best fit of the QSO E_X,Y versusA_V,DL (straight line) in Fig. B.2 agree remarkably well, supports the validity of the preceding analysis.

Once we have computed E_X,Y for all the band pairs and U_min, we remove the QSOs that are considered as outliers in any of the computations to obtain a cleaner sample of good QSOs. We reiterate the full procedure twice using the good QSO sample form the previous iteration, resulting in a final cleanest sample containing 224 245 QSOs with ζ< 3.35 (for which we have Lyα free photometry in the r-, i-, and z-bands), 135,953 with ζ< 2.31 (where we can use the r-band), and 77 633 QSO with ζ< 1.64, where we can use all the SDSS bands. We have an estimate of the intrinsic colours C_X,Y, and an estimate of the reddening E_X,Y for each QSO that is retained by the redshift constraints.

Even though the unobscured QSO intrinsic colours are computed independently for each band pair, we do obtain consistent results across the band pairs, i.e. (B.4)holds for all the bands X, Y, and Z, over all the redshifts ζ considered. Working with the H i column density maps as an estimate of the extinction instead of theA_V,DL gives very similar estimates of C_{X − Y}(ζ), and is independent of any dust modelling, so this means we did not translate potential dust modelling systematics into our QSO estimates.

In order to compare theA_V,DL estimate with a QSO estimate, we need to derive a QSO extinction A_V from the different colour excess E_X,Y. We proceed as follows.

For a given QSO spectrum and extinction curve shape, we can compute the SDSS magnitude increase per dust extinction A_X/A_V for X = u,g,r,i, and z. These ratios depend on the assumed extinction curve and QSO spectral shape, and therefore on the QSO redshift. Using the QSO composite spectrum of Vanden Berk et al. (2001) and the extinction curve presented by Fitzpatrick (1999) parametrized via R_V, we compute the ratios A_X/A_V: (B.5)Figure B.3 shows δX(ζ,R_V = 3.1) as a function of the QSO redshift ζ, for the different bands X = u,g,r,i, and z. Even though the QSO intrinsic colours are strong functions of its redshift, the extinction curves are smooth enough that δX(ζ,R_V = 3.1) is mostly redshift independent.

Using the extinction curves with R_V = 3.1 (which was also used to constrain the optical properties of the grains used in the DL model), for each redshift ζ, we define: (B.6)and (B.7)Finally, for each QSO, we define itsA_V,QSO as the average of the A_{V,QSO, [ X,Y ]} values for all the band pairs that are allowed by its redshift ζ.

Appendix C: Impact of the CIB anisotropies and instrumental noise on the parameter estimation

We study the impact of CIB anisotropies (CIBA) and instrumental (stochastic) noise in our mass estimates in the diffuse ISM (where their effect should be the largest). We simulate data by adding CIBA and instrumental noise to DL SEDs, and fit them with the same technique as we use to fit the observed data. The results quantify the deviations of the recovered parameters from the original ones.

We start by a family of four DL SEDs with U_min = 0.4, 0.6, 0.8, and 1.0, a typical f_PDR = 0.05, and q_PAH = 0.03. We normalize each SED to the mean A_V found for the QSO lines of sight in each U_min. We replicate each SED 100 000 times, add CIB anisotropies and instrumental noise. The noise added has 2 components. We add (band-to-band) independent noise to simulate stochastic instrumental noise with amplitudes given by Pl-MBB, Table B.1, 30′ resolution. We further add a typical CIB SED (also from Pl-MBB, Table B.1, 30′ row), that is completely correlated across the Planck bands, and partially correlated with the IRAS bands, as recommended in Pl-MBB, Appendix B. We finally fit each simulated SED with DL model, as we did in the main data fit.

Figure C.1 shows the recovered Σ_Md divided by the original Σ_Md, and recovered U_min for the SEDs. Each set of points correspond to the different original U_min. The inclined solid line corresponds to the renormalization curve given by Eq. (9), (rescaled to match the mean A_V of the simulated SEDs). There is not a global bias in the recovered Σ_Md, nor U_min; the distribution of the recovered Σ_Md and U_min are centered in the original values. Although CIBA and instrumental noise do generate a trend in the same direction as the renormalization, their impact is significantly smaller than the observed renormalization: they do not span the full range found over the QSOs lines of sight. Moreover, the renormalization found in Sect. 6.3 is independent of the modelling resolution; one obtain similar renormalization coefficients working at 5′, 30′, and 60′ FWHM. For those resolutions, the instrumental noise and CIBA have a very different magnitude, and therefore, their impact would be quite different. Therefore, CIBA and instrumental noise are not a significant source of the A_V systematic departures with respect to U_min found in Sect. 6.3.

Appendix D: Maps in the Planck Legacy Archive

The maps of the dust model parameters, the dust extinction and the model predicted fluxes described in this paper can be obtained from the Planck Legacy Archive (PLA²³). The maps are all at 5′(FWHM) angular resolution in the HEALPix representation with N_side = 2048.

For each quantity, but the χ² of the fit per degree of freedom, there are 2 maps corresponding to the value presented and the corresponding uncertainty. Available maps include our best estimate of the dust extinction (the renormalized A_V,RQ) expressed in magnitude units, and the best fit DL parameters: the dust mass surface density Σ_Md expressed in M_⊙ kpc^-2 units, the starlight intensity heating the bulk of the dust, U_min in units of the ISRF estimated by Mathis et al. (1983) for the solar neighbourhood, the fraction of the dust luminosity from dust heated by intense radiation fields, f_PDR in Eq. (2), which is a dimensionless number between 0 and 1, and the dust mass fraction in small PAH grains q_PAH. We also provide the DL model predicted fluxes in the Planck, IRAS 60 and 100, and WISE 12 bands. Additional information about the file names and the data format is available in the Planck explanatory supplement²⁴.

All Tables

All Figures