© ESO, 2015

1. Introduction

Interstellar dust is an important component of the interstellar medium (ISM): it is studied for its role in the physics and chemistry of the ISM (D’Hendecourt et al. 1985; Mathis 1990; Bakes & Tielens 1994; Le Bourlot et al. 2012), for its role as a tracer of gas column density (Boulanger et al. 1996; Rowles & Froebrich 2009; Planck Collaboration XI 2014) and of magnetic field structure (Chapman et al. 2011; Poidevin et al. 2013; Berdyugin et al. 2014; Planck Collaboration Int. XIX 2015), and for its effect as a foreground in studies of the cosmic microwave background (Planck Collaboration Int. XXII 2015; Planck Collaboration Int. XXX 2015). Many dust models have been made to reproduce the main dust observables, which include extinction curve and albedo, spectral energy distribution (SED) from the near-infrared to the microwave continuum, and elemental abundance constraints (Zubko et al. 2004; Draine & Li 2007; Compiègne et al. 2011; Jones et al. 2013).

Our understanding of interstellar dust, however, is still incomplete. One of the issues faced by the models is that dust is not the same everywhere in the ISM. As revealed by numerous recent observations (Pagani et al. 2010; Köhler et al. 2011; Planck Collaboration XXIV 2011; Planck Collaboration XXV 2011; Köhler et al. 2012; Martin et al. 2012; Paradis et al. 2012; Ysard et al. 2013), dust properties vary from the diffuse ISM to molecular clouds. Variations are especially notable in the far-infrared opacity, i.e. the far-infrared optical depth per unit dust column density τ/N_H. The dust optical depth τ_λ is the function that, in the optically thin limit, modulates the SED: I_λ = τ_λ·B_λ(T), where B_λ(T) is the black-body emission at a temperature T¹.

Variations in dust far-infrared opacity can now be traced in the diffuse ISM itself thanks to the survey of the Planck satellite (Tauber et al. 2010). The far-infrared opacity shows significant variations in the diffuse ISM, although those variations are smaller than the variations observed between the diffuse ISM and molecular clouds (Bot et al. 2009; Planck Collaboration Int. XVII 2014; Planck Collaboration XI 2014).

Planck Collaboration Int. XXIX (2014) examines the dust emission obtained from Planck , IRAS (Neugebauer et al. 1984) and WISE (Wright et al. 2010) data, fitting it over the whole sky with the Draine & Li (2007) dust model, which uses constant optical properties². In this paper, our analysis concentrates on more than 200 000 lines of sight towards quasi-stellar objects (QSOs) observed in the Sloan Digital Sky Survey (SDSS, York et al. 2000) and for which A_V is available³: the A_V from the QSOs is correlated with the value ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}$ deduced from the Planck and IRAS fit. Planck Collaboration Int. XXIX (2014) find that Galactic dust SEDs can be well fitted by the Draine & Li model (Draine & Li 2007) when one uses as free parameters the dust column density – or, equivalently, the V band extinction A_V – and the intensity U_min of the interstellar radiation field (ISRF) in the diffuse ISM. However, even in the diffuse ISM, this model cannot explain the variations in the SED per unit extinction if the ISRF intensity is the only quantity that varies⁴: the Draine & Li (2007) model can reproduce the dust SED, but it systematically overestimates the observed A_V by a factor that depends on U_min and has a mean value of 1.9 in the diffuse ISM.

We aim to extend the analysis by Planck Collaboration Int. XXIX (2014) to the Compiègne et al. (2011) and Jones et al. (2013) dust models and to quantify variations of their far-infrared emission properties in the diffuse ISM. Our analysis will focus on the big grains in thermal equilibrium (radius ~100 nm). The paper is organized as follows. In Sect. 2 we introduce the observational data and the dust models used to fit these data. In Sect. 3 we present two different methodologies for constraining dust parameters using observations, we show how their results differ and we estimate the variation of dust optical properties in the diffuse ISM. In Sect. 4 we discuss the results, identifying a possible shortcoming of most current dust models. Finally, in Sect. 5 we explore the prospects for future dust models.

2. Data and models

2.1. Data

Planck Collaboration Int. XXIX (2014) provides a family of 20 A_V-normalised SEDs, I_λ/A_V, to be used as constraints for future dust models. The quantity A_V is measured on QSO lines-of-sight: each SED is the average over ~10 000 lines of sight with similar U_min. The SEDs have U_min in a range of 0.42 to 0.98, with an average of 0.66 and a standard deviation of 0.14. The QSO lines of sight sample regions of low A_V with a median value of ~0.10, however, the set of SEDs describe the data over a much larger fraction of the sky (≥70%) with A_V up to 1⁵. The wavelengths coverage is from 60 μm to 2.1 mm (between 5000 and 143  GHz).

We use a selection from the Planck Collaboration Int. XXIX (2014) SEDs. Dust models agree in attributing the V band extinction mainly to big grains, which, being in thermal equilibrium, also dominate emission at λ> 100 μm. To compare the same grain type in extinction and in emission, therefore, we use the 100 μm (3000  GHz) IRAS band and the 350, 550, and 850 μm (857, 545 and 353  GHz) Planck bands. We exclude the 60 μm band, where emission is dominated by stochastically-heated smaller grains, because a single band is not sufficient to test small grains properties: a measure of the extinction in the UV is needed to determine the total mass of small particles, and a full mid- and near-infrared emission spectrum is needed to constrain the small particle temperature and size distribution.

Figure 1 shows a subset of the data: the median SED, with U_min = 0.66; the two extreme SEDs, with U_min = 0.98 and U_min = 0.42; the SEDs closest to being 1σ above and below the median, with U_min = 0.80 and 0.52 respectively. The figure shows that the A_V-normalised SEDs are very similar over the Planck wavelength range, and most of the variations are in the IRAS 100 μm band.

Fig. 1 Subset from the 20 AV-normalised dust SEDs in Planck Collaboration Int. XXIX (2014). Each Umin, i.e. each temperature, is shown in a different colour and line style. For clarity of illustration, only five of the 20 SEDs are shown: the warmest and coldest overall SEDs, the warmest and coldest of the SEDs that are within a standard deviation from the average Umin, and the median SED. Symbols indicate the central band wavelengths; error bars indicate the dispersion for each band flux, while instrumental noise is negligible.

The SEDs have both statistical and systematic uncertainties. What is important in statistical uncertainties is their relative value from one band to another, which defines the statistical weight to be applied to each band in the fitting routine. The statistical errors in each band of I_λ/A_V are negligible, since each SED is an average over more than 10 000 observations. The systematic uncertainties in the dust SED come mainly from the photometric calibration: the calibration uncertainty is important at 100 μm (1σ error bar = 13.6%) and at 350 and 550 μm (1σ = 10%); we decided to neglect it at longer wavelengths where it is much smaller () (Miville-Deschênes & Lagache 2005; Planck Collaboration VIII 2014). We estimated the effect of this uncertainty via Monte Carlo simulations. Specifically, for each observed SED, we realized 1000 simulations of random Gaussian-distributed errors for the 100, 350, and 550 μm bands, with the photometric error at 350 and 550 μm set equal since both channels are calibrated on planets (the uncertainties are dominated by the planet model used, not by noise). We performed a fit on all simulations, obtaining 1000 values of each fit parameter for each SED: the mean of these was taken as the fiducial value of the fit parameter and the standard deviation as the uncertainty. The systematic uncertainty in A_V is estimated by Planck Collaboration Int. XXIX (2014) to be ~15%.

2.2. Dust models

We use three dust models with different characteristics: Draine & Li (2007), with the modifications explained in Aniano et al. (2012), Compiègne et al. (2011), and Jones et al. (2013).

Draine & Li (2007, hereafter DL07) is one of the most popular dust models, often used to estimate dust and gas masses in Galactic and extragalactic environments by fitting dust emission (Draine et al. 2007; Aniano et al. 2012). It comprises a size distribution of big silicate grains with optical properties adapted from observations (Draine & Lee 1984; Li & Draine 2001) and carbonaceous grains whose optical properties are assumed to vary continuously from those of graphite for large sizes (≫ 10⁵ C atoms) to those of polycyclic aromatic hydrocarbons (PAHs) for small sizes ( C atoms).

Compiègne et al. (2011, hereafter C11) uses silicate grains with the same properties as DL07, but carbonaceous grains in C11 are divided into two populations: one of PAHs, and one population of amorphous carbon grains, covering a size range from nm-sized to big grains, with optical properties derived from Zubko et al. (1996).

Jones et al. (2013, hereafter J13) is a core-mantle model with optical properties that are derived from laboratory measurements. The J13 model features big grains of amorphous forsterite-like silicates, with iron inclusions and amorphous aromatic carbon mantles; big amorphous carbon grains with aliphatic cores and aromatic mantles; and small aromatic grains (down to sub-nm sizes), which take the place of PAHs. A unique feature of J13 is that it allows for a certain flexibility in its optical properties, which vary according to the band gap of carbonaceous materials and the size of small grains, as well as possible variations in the mantle thickness. We use the standard model aromatic carbon with a band gap of 0.1 eV.

The dust extinction and emission are computed by interpolating an existing library in the case of DL07 (Aniano et al. 2012) and using the DustEM tool⁶ in the case of C11 and J13 (Compiègne et al. 2011). The observational SEDs were fitted using a χ²-minimising procedure. We assumed an interstellar radiation field (ISRF) of the form ${\mathit{u}}_{\mathit{\lambda }}\mathrm{=}{\mathit{G}}_{\mathrm{0}}\mathrm{\times }{\mathit{u}}_{\mathit{\lambda }}^{\mathrm{MMP}\mathrm{83}}$, where ${\mathit{u}}_{\mathit{\lambda }}^{\mathrm{MMP}\mathrm{83}}$ is the mean ISRF in the solar neighbourhood estimated by Mathis et al. (1983) and G₀ is a dimensionless scale factor. The factor G₀ has the same meaning in the model as the U_min introduced in Sect. 1, but we chose to use a different name to more easily distinguish our results and those of Planck Collaboration Int. XXIX (2014). For the rest of the paper, G₀ is our estimate of the ISRF intensity, while U_min is only used to indicate the SEDs and is functionally synonymous with dust temperature.

3. Methodology and results

We analysed the observations using two different methods. In the first method, we fit the SED to obtain the model dust parameters: the radiance, i.e. the integrated SED intensity, the ISRF intensity G₀, and the dust column density, as represented by A_V. In the second method, we estimate the ISRF intensity from the observed A_V and dust radiance, and then use it to predict the model dependent SED.

3.1. SED fitting

Fig. 2 Fitted parameters ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$ and ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$, for the three models. The colour scale goes from blue (cold) to red (warm). The filled symbols correspond to the SEDs shown in Fig. 1. All points have similar relative uncertainty; the error bars for the lowest Umin are shown. The uncertainties are systematic, so they affect each point in the same way; errors on the two axes are strongly anti-correlated (see Sect. 3.1). The horizontal grey lines show ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}\mathrm{=}\mathrm{1}$ and the 15% uncertainty in the AV normalisation (Sect. 2.1). The slanted grey lines are the least-square fits for each model; they were added to more easily tell apart the three models.

The ISRF intensity ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$ is obtained from the fit of the SEDs. The same fit, since the SEDs are normalised by the QSO-derived A_V, returns the ratio of the model extinction to the observed extinction, ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$. The fit results are shown in Fig. 2; each symbol represents a pair of values (${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$, ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$) for a different model and SED. A perfect model would find a value of ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}\mathrm{=}\mathrm{1}$ (thick horizontal grey line). The thinner horizontal grey lines show our 15% systematic uncertainty in A_V (see Sect. 2.1). This uncertainty is identical for all SEDs, and so it affects the significance of the average ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ of the models, but does not affect either the dependence of ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ on G₀ or the differences between models. The error bars on (${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$, ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$) are obtained from Monte-Carlo simulations, as explained in Sect. 2.1. They are strongly anti-correlated, with an average Pearson correlation coefficient of −0.77 (DL07) to −0.84 (C11 and J13). Since these uncertainties are systematic, the plots may be shifted vertically or horizontally, but their shape remains essentially the same.

The J13 model fits the data well, with an average ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}$ which coincides with the expected value. The C11 model overestimates ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}$ by ~25%. The DL07 model shows the largest discrepancy, with an ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}$ that is overestimated by a factor ~1.8, as already pointed out by Planck Collaboration Int. XXIX (2014). All models, however, show a negative correlation between G₀ and ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ whereas ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ should be unity⁷. For comparison, we also fitted the data with a modified black-body. We obtain similar results, which are presented in Appendix A.

3.2. Recovery of G₀ from dust radiance

Fig. 3 Estimation of G0 from the interpolation of ℛ /AV. The horizontal lines are the value of ℛ /AV for the observed SEDs, the curves from the origin are the ${\mathrm{ℛ}}^{\mathrm{model}}\mathrm{\left(}{\mathit{G}}_{\mathrm{0}}\mathrm{\right)}\mathit{/}{\mathit{A}}_{\mathit{V}}^{\mathrm{model}}$. The abscissae of the intersection points are the estimates for G0 for that SED and model. The same 5 SEDs as Fig. 1 are shown, for clarity of illustration, using the same colour and line style scheme.

The SED fitting yields biased results if the optical properties of the model are not optimal, as was shown by Planck Collaboration Int. XXIX (2014) for the DL07 model. It would be useful to find an alternative way to compare models and observations. Therefore, we decided, in addition to fitting the dust SEDs, to recover G₀ using a procedure based on the dust radiance per unit extinction.

The radiance, ℛ, is the total power emitted by dust grains in thermal equilibrium (Planck Collaboration XI 2014). Through the conservation of energy, ℛ corresponds to the power absorbed by grains in thermal equilibrium and can be considered a “heating power” for the dust, which is only determined by its absorption properties, independent of its far-infrared opacity.

The thermal emission of big grains can be reasonably well fitted with a modified black-body, I_λ = B_λ(T)·τ₀·⁽λ/λ₀^{)− β}, with T the observational dust temperature, β its spectral index and τ₀ its optical depth at the reference wavelength λ₀ (here we use 850 μm, equivalent to a frequency ν₀ = 353  GHz). The integration of a modified black-body yields (Planck Collaboration XI 2014) $\mathrm{ℛ}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}\frac{{\mathit{\sigma }}_{\mathrm{s}}}{\mathit{\pi }}{\mathit{T}}^{\mathrm{4}}{\left(\frac{\mathit{kT}}{\mathit{h}{\mathit{\nu }}_{\mathrm{0}}}\right)}^{\mathit{\beta }}\frac{\mathrm{\Gamma }\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\beta }\mathrm{)}\mathit{\zeta }\mathrm{(}\mathrm{4}\mathrm{+}\mathit{\beta }\mathrm{)}}{\mathrm{\Gamma }\mathrm{\left(}\mathrm{4}\mathrm{\right)}\mathit{\zeta }\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{·}$(1)The functions Γ and ζ are the Gamma and Riemann zeta function, respectively, and σ_s the Stefan-Boltzmann constant.

For a given dust model one can compute its A_V-normalised fluxes in the IRAS and Planck bands, fit a modified black-body to the normalised SED (for λ ≥ 100 μm), and substitute the resulting parameters in Eq. (1). The result, ${\mathrm{ℛ}}^{\mathrm{model}}\mathit{/}{\mathit{A}}_{\mathit{V}}^{\mathrm{model}}$, is a linear function of G₀, and can be inverted to find G₀ from the radiance of an A_V-normalised SED. Figure 3 shows ${\mathrm{ℛ}}^{\mathrm{model}}\mathrm{\left(}{\mathit{G}}_{\mathrm{0}}\mathrm{\right)}\mathit{/}{\mathit{A}}_{\mathit{V}}^{\mathrm{model}}$ for the three models. As we mentioned, ℛ only depends on absorption, and so the estimation of ISRF intensity given by this method, ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$, is not affected by variations in dust far-infrared opacity. Since the method of obtaining ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ assumes that the dust optical properties are fixed in extinction, while the method of obtaining ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$ assumes that optical properties are fixed in both extinction and emission, we expect ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ to be much less biased than ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$. This conclusion is supported by the fact that G₀ estimates become much more model dependent when they include optical properties in extinction: ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ varies by ~15% between models, as opposed to the factor ~2 in ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$.

Fig. 4 Data SEDs (top left) and dust model Iλ/ ℛ, normalised by the median data SED (Umin = 0.66) to highlight the spectral shape variations. The model SEDs are computed using ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ as ISRF intensity (Sect. 3.2). The same five SEDs as Fig. 1 are shown, for clarity of illustration, using the same colour and line style scheme. The observations are plotted in grey behind the models to aid comparison.

The values for ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$, shown in Fig. 3, vary by a factor 1.6 between the two extreme SEDs and by 1.3 if we only consider the SEDs within 1σ of the average (which would represent the “typical” G₀ variation rather than the full range). In comparison, in Fig. 2, where the dust optical properties in the model are fixed and the ISRF is the only source of difference between SEDs, ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$ varies by a factor 2.3 between extremes and 1.6 within the “± 1σ” range.

After obtaining ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$, we use it as a fixed parameter in the computation of the SED of the three dust models. The results are shown in Fig. 4. To highlight spectral shape variations, the SEDs have been normalised by the median data SED and then further divided by their ℛ /A_V. The J13 model reproduces the median SED almost perfectly. Predictions based on the C11 model systematically present an excess emission at short wavelengths and low emission at long wavelengths. This situation is even more drastic in the DL07 model. Furthermore, for all models, the normalised SEDs in Fig. 4 spread over a smaller range than the data: the warm (high U_min) SEDs and the cold (low U_min) ones are closer in temperature for the model than they are for the observations.

The differences shown in Fig. 4 between the J13, C11 and DL07 dust models can be attributed to their optical properties in far-infrared emission, as opposed to their optical properties in extinction⁸. We find that the J13 model (Jones et al. 2013) is a better model for the diffuse ISM emission than that presented in Compiègne et al. (2011) and an even better model than the Draine & Li (2007) model. The J13 model uses optical properties for silicates and amorphous carbon based on lab measurements. The Compiègne et al. (2011) and Draine & Li (2007) models use silicates with far-infrared optical properties extrapolated from mid-infrared astronomical observations (Draine & Lee 1984) and empirically adjusted (Li & Draine 2001) to reproduce the FIRAS spectra of Finkbeiner et al. (1999)⁹: these optical properties pertain to silicates that are contaminated by other materials, including possibly a carbonaceous dust component. Furthermore, Draine & Li (2007) is unique among the models in that its carbon dust is assumed to be graphite. Apparently, “astronomical silicate” and graphite are not emissive enough in the far-infrared and submillimeter.

3.3. Variation in dust opacity

The models shown in Fig. 4 assume constant optical properties and are tuned to reproduce the correct ℛ /A_V for each SED, meaning that spectral shape variations in the modelled SEDs are solely driven by the ISRF. It is clear that the full range of spectral shapes exhibited by the observations (Fig. 4, top left) cannot be reproduced by ISRF variations alone. As mentioned in Sect. 3.2, if the dust optical properties are fixed the variations of ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ are only about half of what is needed to reproduce the data. We need dust with variable optical properties to reproduce the observations.

We start by analysing a simplified scenario in which the dust optical properties in extinction are fixed, i.e. A_V is a proxy for the dust column density, and the SED shapes are explained by the variation in far-infrared opacity. Our aim is to quantify this variation, at least to a first approximation. At the end of this section, we see what changes when the assumption of constant extinction is dropped.

Since the SEDs in our work are normalised by A_V, any determination of the far-infrared optical depth from the data actually yields τ_FIR/A_V. In addition, since A_V is a proxy for the dust column density, the variations of τ_FIR/A_V are identical to those of the far-infrared opacity. Figure 2 shows that ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ increases by 40−50% between the warmest and coldest SEDs, independent of the model. This indicates that in the diffuse ISM, τ_FIR/A_V (and therefore opacity) varies by a factor 40−50% between extremes. When considering just the SEDs within 1σ of the average, we obtain an estimate of the “typical” variation: ~20%.

We can compare this estimate to another, obtained with a different technique. We repeat the analysis of Sect. 3.2 on modified C11 and J13 models in which the far-infrared opacity of silicates and carbon has been multiplied by the same factor, independent of grain size, as shown in Fig. 5. This modification is a purely phenomenological artifice, not based on any physical properties of the materials. We only use this modification to estimate the variations in far-infrared opacity that a physical model should attempt to reproduce. We find that C11 would reproduce all the SEDs within 1σ of the average if we allowed its far-infrared opacity to vary between 1.1 and 1.3 times its standard value (between 1.0 and 1.4 to reproduce the full range of observations). The results are similar for the J13 model with scaling factors 0.9 and 1.1 for the 1σ range of SEDs (0.8 and 1.2 for the full range). Thus, we estimate the variations of far-infrared opacity to be ~20% (40−50%) for the typical (full) range of diffuse ISM.

So far, we have assumed that optical properties in extinction remain the same while the far-infrared opacity varies. If we drop the assumption the conclusion must be modified. When the optical properties in extinction vary, A_V is no longer a proxy for the dust column density: our previous conclusions on the variations of τ_FIR/A_V are still valid, but they no longer apply to far-infrared opacity either. Since A_V is proportional to the extinction cross-section in the V band, this means that the variation range found above applies to the ratio of the far-infrared opacity to the extinction cross-section in the V band.

Fig. 5 Effect of variable far-infrared opacity on Iλ/ ℛ for the C11 (left) and J13 (right) models. The same five SEDs as Fig. 1 are shown, except the median SED, using the same colour and line style scheme. The normalisation is the same as Fig. 4. The observations are plotted in grey behind the models; the standard dust models are the solid lines with filled polygonal symbols and the modified dust models are the dashed lines with empty symbols. Numbers show the multiplicative factor of far-infrared opacity.

4. Discussion

While some models are closer to the mean SED, no model is able to reproduce the full extent of observed variations in the SEDs. We recall that in the three models the optical properties of dust are assumed constant. Our results show that the variation of G₀ alone is not sufficient to reproduce the entire span of spectral shapes presented by the observations. Models that endeavour to explain the emission of the diffuse ISM need to include physical processes that change the far-infrared optical depth per unit A_V, or the ratio of the far-infrared opacity to the V band extinction cross-section.

Variations in size distribution can lead to changes in A_V while having little effect on the far-infrared opacity. However, the resulting variations in I_λ/A_V are small given the constraints imposed on G₀ by ℛ /A_V (See Appendix B). A promising route is to consider variations in the grain composition. One process that could alter the grain composition is the accretion of carbon, from the gas phase or from small particles, on big grains. It is readily possible to consider this process within the J13 model because it is tuned to exactly this kind of scenario: its big grains, both silicate and carbonaceous, include an aromatic carbon mantle, the thickness of which can be modified. In the same model, it is also possible to explore the effect of having iron inclusions partially composed of iron sulphide (Köhler et al. 2014). Ysard et al. (2015) shows that the two scenarios, together with variations in the grain size distribution and in the relative abundance of silicate and carbon grains, can explain most of the observed SED variations. Grain shape and structure also play a role in opacity. It is known, for instance, that including porosity in grains increases the far-infrared opacity, while leaving the properties in extinction relatively constant (Ossenkopf 1993; Stepnik et al. 2003; Köhler et al. 2012). Porous grains are usually considered typical of the dense ISM rather than the diffuse ISM, but it has been suggested (e.g. Martin et al. 2012) that they could migrate from the dense environment to the diffuse ISM. Far-infrared opacity also increases when one employs non-spherical grains (e.g. Siebenmorgen et al. 2014), the existence of which is attested by interstellar polarisation. Determining whether a relation between dust evolution and polarisation exists may be an interesting follow-up analysis of the data.

5. Conclusion and perspectives

We confront three dust models (Draine & Li 2007; Compiègne et al. 2011; Jones et al. 2013) to a family of A_V-normalised SEDs obtained from a combination of the Planck, IRAS, and SDSS surveys probing the Galactic diffuse ISM (Planck Collaboration Int. XXIX 2014) .

The widespread method of fitting the SED to obtain the radiation field intensity G₀ and the dust column density (from A_V) gives results that are very model dependent. The average ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}$ of Jones et al. (2013) is close to the observed value, while that of Compiègne et al. (2011) is about 25% too high and that of Draine & Li (2007) too high by a factor ~1.8. Furthermore, all models present a non-physical mass-temperature anti-correlation and yield biased G₀ estimates.

We develop an alternative G₀ estimate that uses the radiance per unit extinction, ℛ /A_V. This estimate is robust with respect to variations in the far-infrared optical properties and gives less biased results. We use these values to compute the SEDs for the models, finding that the variations in observed SEDs cannot be reproduced without varying dust optical properties. We estimate that, to reproduce the typical variations observed, we need to vary the ratio of the far-infrared opacity to the V band extinction cross-section by ~20% in the typical diffuse ISM (up to 40−50% to explain the full range of observations); at the same time, variations of ~30% are needed in the ISRF intensity (up to ~60% for the full range of observations).

Online material

Appendix A: Modified black-body models

Fig. A.1 Scatter plot of ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ and fitted dust temperature for the modified black-body fit. Two alternative estimates of AV are shown (see text), plus an average value. The horizontal grey lines show ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}\mathrm{=}\mathrm{1}$ and the 15% uncertainty in the AV normalisation (Sect. 2.1). All points have similar relative uncertainty; the bars for the lowest Umin are shown. The uncertainties are systematic, so they affect each point in the same way; errors on the two axes are strongly anti-correlated. The Umin colour scheme is the same as Fig. 2; the filled symbols correspond to the SEDs shown in Fig. 1.

Figure A.1 shows two different estimates of A_V as a function of temperature, obtained from a modified black-body fit to the observed SEDs. Planck Collaboration XI (2014) provides two different empirical relations for computing the dust reddening E(B − V) from the emission. One relation considers E(B − V) to be proportional to the 850 μm dust optical depth (E(B − V) /τ₀ = 1.49 × 10⁴), the other relation considers E(B − V) proportional to the radiance (E(B − V) / ℛ = 5.4 × 10⁵). We calculate E(B − V) in both ways, and then we convert the values to A_V using the average diffuse ISM value for R_V = E(B − V) /A_V = 3.1 and compare the results.

The two estimates differ both by their average value and their trend with temperature: the A_V obtained from τ₀ is about 20% too low on average and decreases with temperature like the physical models (Fig. 2); the A_V obtained from ℛ has a good average value but increases with temperature. By way of comparison we also show the geometric average of the two. This average matches the expected value better, despite having no physical justification. Interestingly, the two estimates implicitly make opposite assumptions: the A_V obtained from τ₀ assumes that the dust optical properties are fixed, or at least τ₀/A_V is fixed; and the A_V from ℛ assumes a fixed G₀, or at least a fixed absorbed power per grain¹⁰. This means that cold dust is more emissive than expected from models with fixed dust properties, but less emissive than expected from models where variable optical properties account for all observed variations.

Appendix B: Effects of grain size distribution

The size distribution of dust grains has a strong effect on extinction, but it is not expected to affect the far-infrared opacity, which only depends on the total volume of the grains in the Rayleigh regime. It is natural, therefore, to consider variations in grain size distribution as a way of varying τ_FIR/A_V. Assessing the effect of grain size variation in a physically realistic way is not straightforward: the physical processes that change grain sizes, e.g. shattering, sputtering, accretion, coagulation, also affect its structure and composition; also, the optical properties of the materials themselves may be size-dependent. Simply varying the grain size distribution in a model is therefore not likely to mimic the actual variations in the ISM, but can still provide interesting qualitative insights. In this Appendix we explore the effects of varying the size distribution using the C11 model. While this model does not fit the average I_λ/A_V as well as J13, it is still close enough to be useful for a differential analysis. Its homogeneous grains and constant optical properties allow us to modify the grain size distribution independent of optical properties.

In C11, grains larger than ~10 nm are distributed according to a power law – n(a) ∝ a^α, where a is the grain radius – with an exponential cutoff above ~150 nm. The parameter that mainly controls the size distribution is the exponent of the power law, α, which is −2.8 for carbonaceous grains and –3.4 for silicate grains. We repeat the procedure of Sect. 3.2 varying α by −0.5 and + 0.5 around its standard value. This changes R_V by −0.7 and + 1.0, respectively; by comparison, Fitzpatrick & Massa (2007) give ~0.3 as the typical 1σ dispersion of R_V in the diffuse ISM.

Fig. B.1 Effect of the grain size distribution on Iλ/ ℛ, as predicted by the C11 model. The SEDs shown are the coldest (blue dot-dashed line), the median (green solid line) and the warmest (red dotted line); normalisation is the same as per Fig. 4. The corresponding observations are plotted in grey behind the models. Larger symbols indicate larger average grain size, and smaller symbols indicate smaller average grain size (see text for details).

The results are shown in Fig. B.1. Varying the size distribution has a small impact on the dust SED, and the range of temperatures reproduced is smaller than that observed despite the large, possibly overestimated, span in α. The figure also shows that models with smaller grains are, surprisingly, colder than models with larger grains, i.e. they have lower 100 μm emission and higher long-wavelength emission. This can be explained simply. In our modelling, the radiation field intensity is not fixed, but is rather derived from the observed radiance per unit extinction ℛ /A_V. Changing the size distribution modifies our estimate for G₀, so that models with smaller grains necessitate a weaker radiation field to satisfy those constraints. The decrease in G₀ thus offsets the temperature increase due to size effects, and even reverses it in the case of the C11 model.

The details of the result presented in this appendix are likely to depend on the dust model and parametrization used. Still, varying the grain size distribution without its corresponding change in the dust optical properties is not likely to explain the observed variations of I_λ/A_V. In Ysard et al. (2015) a similar study uses the J13 model, which is instead adapted to reproduce the interplay of grain size and optical properties.

Acknowledgments

F.B. and G.A. acknowledge support from the ERC grant No. 267934.

References

All Figures

Fig. 1 Subset from the 20 AV-normalised dust SEDs in Planck Collaboration Int. XXIX (2014). Each Umin, i.e. each temperature, is shown in a different colour and line style. For clarity of illustration, only five of the 20 SEDs are shown: the warmest and coldest overall SEDs, the warmest and coldest of the SEDs that are within a standard deviation from the average Umin, and the median SED. Symbols indicate the central band wavelengths; error bars indicate the dispersion for each band flux, while instrumental noise is negligible.

In the text

Fig. 2 Fitted parameters ${\mathit{G}}_{\mathrm{0}}^{\mathrm{fit}}$ and ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$, for the three models. The colour scale goes from blue (cold) to red (warm). The filled symbols correspond to the SEDs shown in Fig. 1. All points have similar relative uncertainty; the error bars for the lowest Umin are shown. The uncertainties are systematic, so they affect each point in the same way; errors on the two axes are strongly anti-correlated (see Sect. 3.1). The horizontal grey lines show ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}\mathrm{=}\mathrm{1}$ and the 15% uncertainty in the AV normalisation (Sect. 2.1). The slanted grey lines are the least-square fits for each model; they were added to more easily tell apart the three models.

In the text

Fig. 3 Estimation of G0 from the interpolation of ℛ /AV. The horizontal lines are the value of ℛ /AV for the observed SEDs, the curves from the origin are the ${\mathrm{ℛ}}^{\mathrm{model}}\mathrm{\left(}{\mathit{G}}_{\mathrm{0}}\mathrm{\right)}\mathit{/}{\mathit{A}}_{\mathit{V}}^{\mathrm{model}}$. The abscissae of the intersection points are the estimates for G0 for that SED and model. The same 5 SEDs as Fig. 1 are shown, for clarity of illustration, using the same colour and line style scheme.

In the text

Fig. 4 Data SEDs (top left) and dust model Iλ/ ℛ, normalised by the median data SED (Umin = 0.66) to highlight the spectral shape variations. The model SEDs are computed using ${\mathit{G}}_{\mathrm{0}}^{\mathrm{ℛ}}$ as ISRF intensity (Sect. 3.2). The same five SEDs as Fig. 1 are shown, for clarity of illustration, using the same colour and line style scheme. The observations are plotted in grey behind the models to aid comparison.

In the text

Fig. 5 Effect of variable far-infrared opacity on Iλ/ ℛ for the C11 (left) and J13 (right) models. The same five SEDs as Fig. 1 are shown, except the median SED, using the same colour and line style scheme. The normalisation is the same as Fig. 4. The observations are plotted in grey behind the models; the standard dust models are the solid lines with filled polygonal symbols and the modified dust models are the dashed lines with empty symbols. Numbers show the multiplicative factor of far-infrared opacity.

In the text

Fig. A.1 Scatter plot of ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}$ and fitted dust temperature for the modified black-body fit. Two alternative estimates of AV are shown (see text), plus an average value. The horizontal grey lines show ${\mathit{A}}_{\mathit{V}}^{\mathrm{fit}}\mathit{/}{\mathit{A}}_{\mathit{V}}\mathrm{=}\mathrm{1}$ and the 15% uncertainty in the AV normalisation (Sect. 2.1). All points have similar relative uncertainty; the bars for the lowest Umin are shown. The uncertainties are systematic, so they affect each point in the same way; errors on the two axes are strongly anti-correlated. The Umin colour scheme is the same as Fig. 2; the filled symbols correspond to the SEDs shown in Fig. 1.

In the text

Fig. B.1 Effect of the grain size distribution on Iλ/ ℛ, as predicted by the C11 model. The SEDs shown are the coldest (blue dot-dashed line), the median (green solid line) and the warmest (red dotted line); normalisation is the same as per Fig. 4. The corresponding observations are plotted in grey behind the models. Larger symbols indicate larger average grain size, and smaller symbols indicate smaller average grain size (see text for details).

In the text