© ESO, 2011

1. Introduction

The infrared (IR) spectral energy distribution (SED) is widely used to derive the global properties of a system, such as its instantaneous star formation rate, its dust and eventually gas masses, and the compactness of the star forming region. The advent of the Herschel Space Observatory has opened the most important spectral window to perfect these diagnostics, by observing the far-IR to submillimeter (submm) wavelengths (60–600 μm). Indeed, this regime samples the peak and Rayleigh-Jeans wing of the dust emission. It consequently constrains the emission by grains in thermal equilibrium with the radiation field, present in the different phases of the interstellar medium (ISM), including the coldest, most massive components (down to dust temperatures of T_dust ≳ 12 K). This spectral domain is therefore crucial to derive accurate dust masses, and physical conditions, and can be used as a powerful, unprecedented tool to probe interstellar matter in regions where no other counterpart is accessible.

Unfortunately, there are several fundamental systematic unknowns inherent to dust modelling, which are questioning the reliability of these diagnostics. First, the microscopic properties of the grains are still poorly known. In the Milky Way, the most complete and accurate models are constrained by observations of high latitude cirrus clouds: their IR emission, ultraviolet (UV)-to-near-IR extinction, and elemental depletions (Zubko et al. 2004; Draine & Li 2007; Compiègne et al. 2011). The authors performing these models derive the size distribution and abundance of the different grain species – silicates, carbon grains (graphite or amorphous carbons) and polycyclic aromatic hydrocarbons (PAH). Zubko et al. (2004) demonstrated an important degeneracy by presenting complete fits of the same data set (Galactic emission, extinction and depletion), with five different dust compositions, alternating bare and coated grains, as well as crystalline and amorphous solids. Thus, the derived dust properties depend on the assumed chemical composition of each species. The UV to millimetre (mm) opacities are sensitive to the grain composition. They are derived from sparse constraints including astrophysical features, laboratory spectra of analogs of interstellar dust materials, and theoretical solid state physics (e.g. Weingartner & Draine 2001; Draine 2003b). Their universality is doubtful.

The second major source of uncertainties concerns the macroscopic variations of these microscopic grain properties, as a function of the environment. These variations are numerous; some are speculative:(i) PAHs are known to be destroyed in H ii regions (e.g. Madden et al. 2006); (ii) the variations in the R_V parameter of the extinction curve is interpreted as variations of the grain size distribution (Draine 2003a; Fitzpatrick & Massa 2005); (iii) coagulation occurs in dense regions (e.g. Stepnik et al. 2003; Berné et al. 2007); (iv) blast waves are responsible for grain fragmentation and erosion in the low-velocity phase (Jones et al. 1996) and destruction close to the remnant (Reach et al. 2002); (v) the dust abundances and properties are thought to evolve with the metallicity of the ISM (Galliano et al. 2003, 2005, 2008a); (vi) a transition from amorphous to crystalline silicates is observed in protostellar objects (e.g. Hallenbeck et al. 2000; Poteet et al. 2011). This list is not exhaustive. Moreover, when considering the SED of a given region, it is possible to confuse variations of the physical properties of the grains (e.g. their optical properties) with variations of their physical conditions (e.g. the starlight intensity to which they are exposed). This problem becomes even more intricate, when considering the integrated SED of a galaxy.

The various processes controlling the lifecycle of dust throughout the ISM are not known with enough precision to break these kinds of degeneracies. Even the origin of interstellar dust is uncertain. The contribution of supernovae (SN) and asymptotic giant branch (AGB) stars to the observed content of ISM dust, and the dust growth in interstellar cloud is still debated (e.g. Galliano et al. 2008a; Draine 2009). Dust is believed to constantly evolve throughout the ISM, being photoprocessed, altered by cosmic rays, accreting atoms in dense regions, and being shattered in shocks (Jones 2004, for a review). We are compelled to find observational cases where there will be enough redundancy in the data to isolate one of these processes. This is the goal of this paper.

The present article scrutinizes the different methodological biases, as well as the fundamental physical processes affecting the dust mass estimate in galaxies. Our demonstration is performed on the Herschel and Spitzer observations of the Large Magellanic Cloud (LMC; d = 50 ± 2.5 kpc; Schaefer 2008). Due to its proximity, it is an ideal laboratory to study the variations of the far-IR properties, down to spatial scales of  ≃10 pc (SPIRE_500 μm angular resolution of 36′′). Moreover, it offers an environment containing massive star clusters, allowing us to study the impact of intense star formation on the surrounding ISM. Finally, its metallicity is moderately sub-solar, with (O/H)_LMC ≃ 0.5 × (O/H)_⊙ and (C/H)_LMC ≃ 0.3 × (C/H)_⊙ (Pagel 2003). The comparison of its dust properties with those of the Galaxy therefore provides insights on cosmic dust evolution.

In general, low-metallicity dwarf galaxies, a category to which the LMC belongs, have peculiar dust properties. They exhibit a deficit of PAH strength, that appears to be correlated with the metallicity of their ISM. The origin of this trend is still debated: (i) PAHs could be more massively destroyed by permeating hard radiation, in sub-solar ISM (e.g. Madden et al. 2006); (ii) the delayed injection of PAHs by AGB stars could explain their lower intrinsic abundance, in young systems (Galliano et al. 2008a); (iii) the PAHs could form in molecular clouds, which have a lower filling factor, in low-metallicity environments (e.g. Sandstrom et al. 2010).

The IR SED of dwarf galaxies usually peaks at shorter wavelengths, indicating hotter equilibrium grains, on average. In addition, the mid-IR continuum is steeply rising, similarly to what is observed in Galactic compact H ii regions (e.g. Peeters et al. 2002). This peculiar mid-IR continuum, was modelled by Galliano et al. (2003, 2005) with an increase of the very small grain abundances, which could be the consequence of the high number density of shock waves. It is consistent with the peculiar shape of the extinction curve, in the Magellanic clouds: their lower 2175 Å bump, and their steeper near-UV rise could be the result of an excess of small grains (Weingartner & Draine 2001; Galliano et al. 2003, 2005). This typical continumm shape was reported as a “mid-IR excess” compared to the Galactic SED by Bernard et al. (2008), in the LMC.

The dust-to-gas mass ratio increases with metallicity (Lisenfeld & Ferrara 1998; Draine et al. 2007; Galliano et al. 2008a; Engelbracht et al. 2008). At first order, it can be understood since dust is made out of the metals synthesized by the various stellar populations. However, the detailed dependency of the dust-to-gas mass ratio with the metallicity remains unknown. In addition, the gas mass, especially the molecular phase, is uncertain, because the CO-to-H₂ conversion factor is a strong unknown function of the environment – in particular, it is a function of the metallicity. A given CO line intensity will translate in a larger molecular gas mass in a lower metallicity environments (e.g. Madden et al. 1997; Israel 1997; Leroy et al. 2007, 2009, 2011). Moreover, we can expect the presence of large quantities of H₂ in regions where no CO at all is detected (Madden 2000). Bernard et al. (2008) unveiled a “far-IR excess”, compared to the gas column density, in the LMC, which is likely the evidence of such a gas reservoir.

Finally, at submm wavelengths, the SED of dwarf galaxies differs significantly from solar metallicity systems. Galliano et al. (2003, 2005) reported in four blue compact dwarf galaxies an excess emission at 850 μm (SCUBA) and 1.2 mm (MAMBO). Such an excess was then confirmed in other similar systems (Dumke et al. 2004; Galametz et al. 2009, 2010; Bot et al. 2010; Grossi et al. 2010). It extends up to centimetre (cm) wavelengths in the LMC and SMC (Israel et al. 2010; Planck Collaboration et al. 2011a). Although, the COBE data of our Galaxy presented a submm excess (Reach et al. 1995), the intensity of this excess is much more pronounced in low-metallicity systems. Several explanations are in competition, for the origin of this excess: (i) very cold dust, in dense clumps, accounting for a large fraction of the dust budget of the galaxy (Galliano et al. 2003, 2005; Galametz et al. 2009; O’Halloran et al. 2010); (ii) temperature dependent submm emissivity (Meny et al. 2007); (iii) rapidly spinning grains in addition to another component (Bot et al. 2010; Planck Collaboration et al. 2011a).

The first Herschel observations of the LMC showed that the slope of the submm SED was flatter than in the Galaxy (Gordon et al. 2010). Meixner et al. (2010) showed that modelling this SED with standard Galactic grain properties required too much mass, and therefore concluded that it required modified grain optical properties. Roman-Duval et al. (2010) confirmed the Bernard et al. (2008) “far-IR excess” toward several molecular clouds. Hony et al. (2010) demonstrated the complex structure of two massive starforming regions. Gordon et al. (2010) reported a SPIRE_500 μm excess, which is likely the rise of the submm excess previously discussed.

The unprecedented sensitivity and wavelength coverage of Herschel, at far-IR/submm wavelengths, allow us, for the first time, to study in detail processes that were previously glimpsed at. With a rigorous method, accounting for the different sources of error, it is now possible to unveil the systematic effects inherent to SED modelling. The common assumptions, concerning the universality of dust properties, the accuracy of gas mass estimates, and the homogeneity of gas-to-dust mass ratios can be confronted by data. From a technical point of view, all the processes that have been previously described here define the required model parameters, as well as the unknowns when interpreting the IR/submm emission of the LMC.

For that purpose, the present paper is organized as follows. Section 2 presents the data set upon which we have based our analysis, and discusses the reference observational quantities we consider. In Sect. 3, we present our SED model, and the way we rigorously propagate the various sources of observational errors through the entire fitting procedure, in order to provide reliable errors on the derived physical parameters. Section 4 attempts to reconcile different interpretations of the peculiar far-IR properties of the LMC: modified grain composition and/or undetected gas reservoir. We end by a discussion on the origin of the SPIRE_500 μm excess. Section 5 synthesizes the paper and emphasizes the consequences of our findings. The various appendices give details on technical points that would otherwise alter the flow of the discussion, if they were included in the main text.

2. The data set

The data set we are using are the science demonstration (SD) Herschel/SPIRE observations of the LMC (Meixner et al. 2010), together with its Spitzer/IRAC and Spitzer/MIPS data (Meixner et al. 2006), and some ancillary data. These SD Herschel data cover only one fourth of the LMC. Although the complete PACS and SPIRE maps have now been obtained, we have performed our analysis on the sole SD strip, since we want to demonstrate general effects on the dust estimate, that do not require the totality of the LMC.

2.1. Herschel observations and data reduction

We use the SPIRE maps presented by Meixner et al. (2010). We point the reader to this paper for the detailed description of the data reduction. The SD SPIRE data of the LMC cover a 2° × 8°, at 250, 350 and 500 μm. The extended source calibration was performed assuming SPIRE beam areas of 395, 740 and 1517″², at 250, 350 and 500 μm, respectively.

A background was subtracted, taking as a reference the two outer edges of the strip. Those edges are supposed to be out of the LMC. The same regions are considered for the background subtraction at other wavelengths, and for the gas maps.

First, as discussed by Bernard et al. (2008), there are residual foreground Galactic filamentary structures. We used the H i map, whose velocity range corresponds to the Galaxy, in order to quantify the contribution of these fluctuations. Converting the Galactic H i column density into Galactic IR emission (using the Zubko et al. 2004, model), we find that the (non-subtracted) foreground accounts for  ≃15% of the IR power of the strip. When this foreground is subtracted, the remaining fluctuations are on average  ≃1% of the IR power. Therefore, this contamination is smaller than our uncertainties.

Second, we note that the larger scale Planck images (Planck Collaboration et al. 2011a) show that the edges of our maps are not completely beyond the LMC. The Planck images show that there is outer emission, which is colder than that of the rest of the LMC, in particular to the South of the strip. This oversubtraction may bias the dust temperatures, that we will derive in Sect. 4, toward hotter emission. However, as will be discussed in Sect. 4, our results are based on an excess of cold emission in the SPIRE bands. Therefore, this potential oversubtraction is conservative.

We have compared the data of the SD paper, which were a single scan, with the more recent, accurately calibrated data, which also include a cross scan. The relative difference between the integrated SPIRE fluxes of the two sets is less than 2%.

2.2. Spitzer data

The Spitzer data include the four IRAC bands (3.6, 4.5, 5.8, 8.0 μm) and the three MIPS bands (24, 70 and 160 μm). They have been presented by Meixner et al. (2006).

2.3. Gas tracers

To compare the gas and dust contents in the LMC, we have completed our data set with observations of the atomic and molecular gas phases. These maps were presented by Roman-Duval et al. (2010).

We use the [H i]_21 cm map observed by Kim et al. (2003). Their original beam size is 1.0′. The total atomic gas mass in the strip is ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{=}\mathrm{7.14}\hspace{0.17em}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{M}}_{\mathrm{\odot }}$. The noise at the original resolution is σ_H   i(1.0′) ≃ 1.07 M_⊙   pc^-2.

We use the ¹²CO(J = 1 → 0)_2.6   mm map observed by Fukui et al. (2008, with the NANTEN telescope). The beam size is 2.6′. We assume a constant CO-to-H₂ conversion factor of X_CO = 7 × 10²⁰H   cm^-2   (K   km   s^-1)^-1, based on the virial estimate of Fukui et al. (2008). It accounts for the variation of the X_CO factor with metallicity, due to less efficient shielding of CO by dust and to a lower intrinsic C and O abundance We emphasize that this value of X_CO is larger than in the Milky Way. It accounts for discrepancies of the CO-to-H₂ conversion factor in regions where CO is detected. However, it does not take into account potential regions where large envelopes of H₂ could be present, but where the CO would be massively photodissociated, and therefore not detected in emission by ground based radio telescopes. With this conversion, the total molecular gas mass is ${\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\mathrm{2.23}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{M}}_{\mathrm{\odot }}$. The noise at the original resolution is 6.05 M_⊙   pc^-2. The total uncertainty in ${\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}$ is dominated by the uncertainty in X_CO itself. We will discuss that point in Sect. 4.3.

The gas masses above include the mass of helium and heavier elements. The mean atomic weight used is: ${\mathit{\mu }}_{\mathrm{LMC}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}{\mathit{Y}}_{\mathrm{\odot }}\mathrm{-}{\mathit{Z}}_{\mathrm{LMC}}}\mathrm{=}\mathrm{1.34}\mathit{,}$(1)where the mass fractions of helium and heavy elements are Y_⊙ = 0.248 and Z_LMC = 0.5 × Z_⊙ = 8.5 × 10^-3, respectively (Grevesse & Sauval 1998; Pagel 2003). The fraction of molecular gas in the strip is ${\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathit{/}\mathrm{(}{\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{)}\mathrm{\simeq }\mathrm{24}\mathrm{\%}$. It is higher than integrated over the entire LMC (≃10%; Bernard et al. 2008), since the strip includes a large number of molecular clouds.

2.4. Exploring the effects of spatial resolution

All our maps are regridded and reprojected on a common frame. They have been convolved with various kernels, in order to match the spatial resolution of MIPS_160 μm (38″). The full process is described in Gordon et al. (2010).

We aim at studying the systematic effects that would bias the dust mass estimates. The spatial resolution is one of these effects. Indeed, SED models are highly non-linear, since the power emitted by a grain at thermal equilibrium with the radiation field (temperature T_dust) is proportional to ${\mathit{T}}_{\mathrm{dust}}^{\mathrm{4}\mathrm{+}\mathit{\beta }}\mathrm{\simeq }{\mathit{T}}_{\mathrm{dust}}^{\mathrm{6}}$, where β ≃ 2 is the standard “emissivity index” (Eq. (A.1)). Therefore, the sum of the modelled dust masses of N regions is likely to be different than the modelled dust mass of the sum of the emissions of these N regions. In order to study this effect, we will model several maps of the same region, but with different pixel sizes.

Table 1 lists the different resolutions. The highest spatial resolution (R1, 42″, 10 pc) is slightly larger than the largest beam size (MIPS_160 μm). We then construct each set of maps (Spitzer, Herschel and gas) by summing the pixels in a 2 × 2 pixel window. We repeat this process until we reach the size of the full integrated strip (R10). The resolution of the combined gas maps is lower than the dust maps. We will therefore not study the spatial distribution of gas-to-dust mass ratios at resolutions higher than R4. From a technical point of view, the ¹²CO(J = 1 → 0)_2.6   mm map does not cover the entire strip. Consequently, we define a subsample of pixels where both the gas and dust maps are defined.

3. A phenomenological dust SED model

3.1. Motivations

We have developed a model aimed at accurately fitting the observed mid-IR to mm SEDs of various regions of the ISM of the LMC. At the spatial resolutions we consider here (≳ 10 pc), the SEDs will likely be the combination of several regions with different physical conditions – photodissociation regions (PDR), diffuse ISM, etc. In principle, we should perform a radiative transfer model, in order to determine the irradiation of each mass element of the ISM, then compute its spectrum, and transfer the IR radiation to the observer. However, we do not have the necessary information on the detailed matter and stellar 3D distributions, at these spatial scales, to constrain this type of model. Moreover, such an analysis is unnecessary in our case. Indeed, we are interested in the dust mass estimate. The mass is dominated by large grains, at thermal equilibrium with the radiation field. The spectrum of these grains depends only on the stellar power they absorb (or on their equilibrium temperature), and not on the details of the stellar spectrum and spatial distribution.

This is not going to be true for grains which are out of thermal equilibrium with the radiation field (with typical radius a ≲ 10 nm), especially PAHs. The spectrum of these grains depends on the hardness of the radiation field which determines the maximum temperature up to which the grain is fluctuating. The fact is that most of the emission of these grains arises at short wavelengths (λ ≲ 50μm), and is not contaminating the far-IR-to-submm SED.

Finally, the regions we are considering here are optically thin in the IR. Some compact sources show signs of absorption in the mid-IR (9.8 and 18 μm silicate features, 15.2 μm CO₂ ice feature, etc.; Kemper et al. 2010; Hony et al. 2011, in prep.), but the bulk of the grain emission, at longer wavelengths is unaffected.

Considering the previously exposed arguments, we could derive reliable dust masses by simply fitting the observed SEDs with a combination of several modified black bodies. Nonetheless, we still choose to fit a combination of realistic dust models, even if the very small grain (VSG) and PAH spectra will not be perfectly accurate, due to the lack of constraints on the radiation field hardness. This approach is not providing significantly better mass estimates, but is providing more reliable estimates of the physical conditions of the hottest equilibrium grains. This is crucial for the interpretation, as will be demonstrated in Sects. 4.3 and 4.4. In particular, our approach allows us to avoid unphysical fits where a hot equilibrium component will be fit in place of the PAH emission.

Our phenomenological dust SED model can be decomposed in two levels:

• 1.

  the dust SED of a mass element of the ISM, which is controlledby the microscopic grain properties;

• 2.

  the synthesis of several mass elements, to account for the macroscopic variations of the illumination conditions.

Sections 3.2 and 3.3 details these two levels. This model was previously used, in particular, by Galametz et al. (2009, 2010), O’Halloran et al. (2010), Cormier et al. (2010), Hony et al. (2010) and Meixner et al. (2010).

3.2. Dust SED of a mass element of the ISM

Let’s consider the SED emitted by an element of mass of the ISM, where we can assume that the starlight intensity is uniform. For simplicity, we assume that the starlight intensity heating the grains has the spectral shape of the interstellar radiation field (ISRF) of the solar neighborhood (noted ${\mathit{U}}_{\mathit{\lambda }}^{\mathrm{\odot }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}$; Mathis et al. 1983). We parametrize its integrated intensity by: $\mathit{U}\mathrm{=}\frac{{\mathrm{\int }}_{\mathrm{0.0912}\mathit{\mu }\mathrm{m}}^{\mathrm{8}\mathit{\mu }\mathrm{m}}{\mathit{U}}_{\mathit{\lambda }}^{\mathrm{\odot }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{d}\mathit{\lambda }}{\mathrm{2.2}\mathrm{\times }{\mathrm{10}}^{-5}\mathrm{W}\hspace{0.17em}{\mathrm{m}}^{-2}}\mathrm{·}$(2)The value U = 1 corresponds to the intensity of the solar neighborhood. In these conditions, large interstellar silicates have an equilibrium temperature of  ≃17.5 K. It is possible that the ISRF of the LMC differs from ${\mathit{U}}_{\mathit{\lambda }}^{\mathrm{\odot }}$, since it has younger stellar populations. This spectral shape is also likely to vary spatially within the LMC, being harder in transparent regions, and redder in dense regions. However, as explained previously, this will impact only the PAH and VSG spectra, which are not determinant to estimate the total dust mass.

In the same way, we are adopting the Galactic framework, by using the grain properties of Zubko et al. (2004). We choose the bare grain model, with solar abundance constraints (BARE-GR-S). The abundance and size distribution of each grain species was constrained by fitting the IR emission, UV-visible extinction and elemental depletions of the Galactic diffuse ISM. We have updated this model with the new Draine & Li (2007) PAH optical properties, that includes more accurate band profiles, based on Spitzer spectra.

The optical properties and enthalpies considered here are summarized in Table 2. The grain cross-sections are computed using a Mie code, and following the method of Laor & Draine (1993, Sect. 2.2; Appendix A. The temperature fluctuations are computed for each grain size of each component, and for each starlight intensity, using the transition matrix method (Guhathakurta & Draine 1989).

The specific monochromatic luminosity of a mass element of ISM exposed to the starlight intensity U is: ${\mathit{l}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{PAH}}{\mathit{l}}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}\mathrm{+}{\mathit{f}}_{\mathrm{carb}}{\mathit{l}}_{\mathit{\nu }}^{\mathrm{carb}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}\mathrm{+}{\mathit{f}}_{\mathrm{sil}}{\mathit{l}}_{\mathit{\nu }}^{\mathrm{sil}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}\mathit{,}$(3)where f_PAH, f_carb and f_sil are the mass fractions of PAHs (charge fraction of 1/2), graphite and silicate (f_PAH + f_carb + f_sil = 1), and ${\mathit{l}}_{\mathit{\nu }}^{\mathrm{PAH}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}$, ${\mathit{l}}_{\mathit{\nu }}^{\mathrm{carb}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}$, ${\mathit{l}}_{\mathit{\nu }}^{\mathrm{sil}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}$ are the corresponding size distribution integrated specific monochromatic luminosities. In this paper, we keep the mass fractions to the Galactic values (Zubko et al. 2004), except f_PAH that we vary to fit the observed IRAC_8 μm band.

Throughout this paper, we will systematically compare the two following dust compositions, in order to study a possible evolution of composition between the Milky Way and the LMC.

• 1.

  “The standard model” (hereafter labeled Std) is the originalZubko et al. (2004) Galactic graincomposition, made of PAHs, graphite and silicates. The effectivesubmillimeter opacity of this dust corresponds to an emissivityindex β ≃ 2 (Appendix A).

• 2.

  “The AC model” (hereafter labeled AC) is the “standard model”, but replacing graphite by amorphous carbons (ACAR; Zubko et al. 1996). This substitution is arbitrary. The purpose of this model is to test realistic compositions having a higher effective submillimeter opacity (β ≃ 1.7; Appendix A), without violating the elemental abundances.

We note that the heat capacities of amorphous carbons are unknown. We adopt those of graphite in replacement. It might be a crude approximation. However, this inconsistency will affect only the stochastically heated grains, which do not contribute significantly to the dust mass.

The individual SEDs are shown in Fig. 1. The two models show similar features.

• 1.

  The far-IR peak is dominated by grains in thermal equilibriumwith the radiation field. Their spectrum is roughly a modifiedblack body. The peak wavelength of the emission shifts to shorterwavelengths when U rises, as the equilibrium temperature rises.

• 2.

  The mid-IR continuum, for U ≲ 10⁴, is dominated by small grains (radius a ≲ 0.01μm) and PAHs (prominent emission bands at 3.3, 6.2, 7.7, 8.6, 11.3 μm), both out of equilibrium with the radiation field. These grains are being heated by single photon events and the spectral shape of their emission is independent of U. Their spectrum normalized to L_IR (or U) is constant. The change of shape of the mid-IR spectrum with U is only due to the contribution of equilibrium grains at these wavelengths, when their temperature reaches T_eq ≳ 80 K (U ≳ 10⁴). In particular, the prominent 9.7μm silicate feature in emission dominates the mid-IR wavelengths, in this temperature regime.

3.3. Synthetic multi-environment SED

Each SED we model in this paper is likely to be the combination of the emission from regions with different physical conditions. To account for this diversity of conditions, we make the following assumptions.

• 1.

  We assume that the dust properties are uniform within themodelled region: the size distribution, and mass fractions areconstant. Only the starlight intensity varies. This is anapproximation. We will discuss in Sect. 4.3.6potential local variations of the grain properties.

• 2.

  The distribution of starlight intensities per unit dust mass, throughout the region, can be approximated by a power-law (Dale et al. 2001): $\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{\propto }{\mathit{U}}^{\mathrm{-}\mathit{\alpha }}\mathrm{with}{\mathit{U}}_{\min }\mathrm{\le }\mathit{U}\mathrm{\le }{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\mathit{.}$(4)This is an empirical prescription. Dale et al. (2001, Sect. 5.5) provide a physical justification of this formulation. However, its main advantage is that it allows for flexible parametrizing of the physical conditions. A more complex formulation is discussed in Appendix C.1.

The total dust mass of each modelled region is therefore: ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}{\mathrm{\int }}_{{\mathit{U}}_{\min }}^{{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}}\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{d}\mathit{U}\mathit{.}$(5)In addition, to subtract the stellar contribution from the mid-IR bands, we add a stellar continuum, parametrized by the stellar mass in the region M _⋆ : ${\mathit{L}}_{\mathit{\nu }}^{\mathit{\star }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathit{\star }}\mathrm{\times }{\mathit{l}}_{\mathit{\nu }}^{\mathit{\star }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{,}$(6)where ${\mathit{l}}_{\mathit{\nu }}^{\mathit{\star }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}$ is the specific monochromatic luminosity of a 1 Gyr stellar population, synthesized with the model PEGASE (Fioc & Rocca-Volmerange 1997). Since, this population is constrained mainly by the IRAC_3.6 μm and IRAC_4.5 μm bands, the age of the populations do not have a big effect. On the other hand, this stellar mass is poorly determined and should not be trusted. The only purpose of this component is to give a better χ². To accurately determine the mass of the stellar populations, we would need to take into account shorter wavelengths. This is not the purpose of this paper.

In summary, the total monochromatic luminosity of the model is: $\begin{array}{ccc}{\mathit{L}}_{\mathit{\nu }}^{\mathrm{mod}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\\ {\mathit{U}}_{\min }\end{array}{\mathit{l}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{\left(}\mathit{U,\lambda }\mathrm{\right)}\mathrm{\times }\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{d}\mathit{U}\\ & & \end{array}$(7)For each waveband λ_i, where the observed monochromatic luminosity is ${\mathit{L}}_{\mathit{\nu }}^{\mathrm{obs}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}$, we compute the synthetic photometry, ${\mathit{L}}_{\mathit{\nu }}^{\mathrm{mod}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}$, by convolving the model with the instrumental spectral response, using the appropriate conventions provided by the user’s manuals of each instrument. We minimize the χ², using a Levenberg-Marquart algorithm (Press et al. 1992). The χ² is weighted as follows: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{n}}\frac{{\mathrm{\left(}{\mathit{L}}_{\mathit{\nu }}^{\mathrm{obs}}\mathrm{\right(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{)}\mathrm{-}{\mathit{L}}_{\mathit{\nu }}^{\mathrm{mod}}\mathrm{(}{\mathit{\lambda }}_{\mathit{i}}{\mathrm{\left)}}^{\mathrm{\right)}}}^{\mathrm{2}}}{{\left(\mathrm{\Delta }{\mathit{L}}_{\mathit{\nu }}^{\mathrm{rms}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}\right)}^{\mathrm{2}}\mathrm{+}{\mathrm{\left(}\mathrm{\Delta }{\mathit{L}}_{\mathit{\nu }}^{\mathrm{cal}}\mathrm{\right(}{\mathit{\lambda }}_{\mathit{i}}{\mathrm{\left)}}^{\mathrm{\right)}}}^{\mathrm{2}}}\mathit{,}$(8)where $\mathrm{\Delta }{\mathit{L}}_{\mathit{\nu }}^{\mathrm{rms}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}$ and $\mathrm{\Delta }{\mathit{L}}_{\mathit{\nu }}^{\mathrm{cal}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{i}}\mathrm{\right)}$ are respectively the rms and calibration errors of the waveband centered at wavelength λ_i (see Sect. 3.4). We do not use the SPIRE_500 μm flux as a constraint, because of its excess relative to the model, as discussed by Gordon et al. (2010). Instead, we will study the behaviour of this excess in Sect. 4.4, in order to attempt to decipher its origin.

Figure 2 demonstrates this model on the integrated strip SED (R10). The two panels highlight the degeneracy between the starlight intensity distribution and the submillimeter opacities, by comparing the fit of the two compositions to the same observations. Having a flatter submillimeter opacity gives a fit with less massive, hotter dust populations.

Table 3 summarizes the parameters of the model. In particular, the last column of Table 3 describes the behaviour of the free parameters. When an SED fit is performed, a slight variation of one parameter value, due to observational errors, may systematically be compensated by the variation of another parameter. Although the free parameters are rigorously independent, this effect may induce a correlation between these parameters. The Monte-Carlo error analysis that will be discussed in Sect. 3.4.2 is a good way to quantify these correlations. The most striking example, in our case, is the correlation between the three parameters controlling the starlight intensity distribution (α, U_min and ΔU; Eq. (4)). Moreover, these parameters do not have a physical meaning. Rather than discussing their values, it is more convenient to consider the first two moments of the starlight intensity distribution: $\begin{array}{ccc}\mathrm{⟨}\mathit{U}\mathrm{⟩}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{dust}}}{\mathrm{\int }}_{{\mathit{U}}_{\min }}^{{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}}\mathit{U}\mathrm{\times }\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{d}\mathit{U}\\ {\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{U}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{M}}_{\mathrm{dust}}}{\mathrm{\int }}_{{\mathit{U}}_{\min }}^{{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}}{\left(\mathit{U}\mathrm{-}\mathrm{⟨}\mathit{U}\mathrm{⟩}\right)}^{\mathrm{2}}\mathrm{\times }\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}}{\mathrm{d}\mathit{U}}\mathrm{d}\mathit{U,}\end{array}$which develop into: $\begin{array}{ccc}\mathrm{⟨}\mathit{U}\mathrm{⟩}\mathrm{=}\{\begin{array}{c}\\ & \\ & \\ & \\ & \end{array}& & \end{array}$(11)and: ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{U}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \frac{\mathrm{1}\mathrm{-}\mathit{\alpha }}{\mathrm{3}\mathrm{-}\mathit{\alpha }}\frac{{\left({\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\right)}^{\mathrm{3}\mathrm{-}\mathit{\alpha }}\mathrm{-}{\mathit{U}}_{\min }^{\mathrm{3}\mathrm{-}\mathit{\alpha }}}{{\left({\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\right)}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathrm{-}{\mathit{U}}_{\min }^{\mathrm{1}\mathrm{-}\mathit{\alpha }}}\mathrm{-}\mathrm{⟨}\mathit{U}{\mathrm{⟩}}^{\mathrm{2}}& \mathrm{if}\mathit{\alpha }\ne \mathrm{1}\\ & \\ \frac{\mathrm{1}}{\mathrm{2}}\frac{{\left({\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\right)}^{\mathrm{2}}\mathrm{-}{\mathit{U}}_{\min }^{\mathrm{2}}}{\ln \left({\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}\right)\mathrm{-}\ln {\mathit{U}}_{\min }}\mathrm{-}\mathrm{⟨}\mathit{U}{\mathrm{⟩}}^{\mathrm{2}}& \mathrm{if}\mathit{\alpha }\mathrm{=}\mathrm{1.}\end{array}$(12)We define the infrared luminosity as the power emitted by the dust: ${\mathit{L}}_{\mathrm{IR}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{L}}_{\mathit{\nu }}^{\mathrm{dust}}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{d}\mathit{\nu }\mathit{.}$(13)Finally, the “gas-to-dust mass ratio” is the ratio of the total mass of gas (H i and H₂, with helium and heavy elements) to the total dust mass, in the same region: ${\mathit{G}}_{\mathrm{dust}}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{M}}_{\mathrm{dust}}}\mathrm{·}$(14)For comparison, the Galactic value is (Zubko et al. 2004): ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{\odot }}\mathrm{\simeq }\mathrm{158.}$(15)It is important to note that, for the discussion in Sect. 4, this value is consistent with the dust properties we use. Moreover, this model is consistent with the elemental depletion constraints.

3.4. Rigorous error propagation

Since our dust model is highly non-linear, it is crucial to rigorously propagate the observational errors through the entire fitting procedure, taking into account the fact that some errors are independent and others are correlated. In this way, we will be able to quote consistent errors on the parameters. We first need to identify the various sources of error

3.4.1. Sources of observational error

For each wavelength, at each spatial resolution, we measure the noise of the map by taking the standard deviation of the pixel values in what is considered to be the background, i.e. the upper and lower ends of the strip. These values are given in Table 4. These errors are independent from one wavelength to the other, and from one pixel to the other. A simple check of the distribution of pixel values shows that the uncertainty is well described by a Gaussian.

The calibration error is the error on the flux conversion factor. This error is therefore correlated between each pixel. It can be synthesized as follows.

• IRAC:

  the 1σ calibration uncertainty is σ_cal(IRAC) = 2% (Reach et al. 2005). The different wavelengths are correlated.

• MIPS_24 μm:

  Engelbracht et al. (2007) quote a 1σ calibration error of σ_cal(MIPS_24   μm) = 4%. It is independent of the other wavebands.

• MIPS_70 μm:

  Gordon et al. (2007) quote a 1σ calibration error of σ_cal(MIPS_70   μm) = 5% for the coarse scale mapping used for the SAGE observations.

• MIPS_160 μm:

  Stansberry et al. (2007) report a 1σ calibration error of σ_cal(MIPS_160   μm) = 12%. This error is correlated with the MIPS_70 μm error.

• Although Swinyard et al. (2010) report a calibration error of σ_cal(SPIRE) = 15%, it is necessary to decompose this error into its components (SPIRE consortium 2010), as the SPIRE fluxes are the most crucial constraints on the dust mass and emissivity.

  • 1.

    The 3σ error on the calibration model is $\mathrm{3}\mathrm{\times }{\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{model}}\mathrm{\left(}\mathrm{SPIRE}\mathrm{\right)}\mathrm{=}\mathrm{5}\mathrm{\%}$, each waveband being correlated.

  • 2.

    The noise in the calibration observations (Ceres) are:

    • ${\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{noise}}$(SPIRE_250 μm) = 7%;

    • ${\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{noise}}$(SPIRE_350 μm) = 12%;

    • ${\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{noise}}$(SPIRE_500 μm) = 6%.

    These errors are independent.

  • 3.

    The error on the beam area are ${\mathit{\sigma }}_{\mathrm{cal}}^{\mathrm{beam}}\mathrm{\left(}\mathrm{SPIRE}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}\mathrm{\%}\mathit{/}\mathrm{+}\mathrm{5}\mathrm{\%}$ and they are independent.

The latter SPIRE calibration errors correspond to the initial calibration performed during the SD phase. We kept this error, in order to remain consistent with our data set. In any case, this estimate is conservative, the new calibration errors being smaller.

3.4.2. Monte-Carlo iterations

We propagate the observational errors detailed in Sect. 3.4.1, by performing Monte-Carlo iterations of each fit. More precisely, for each observed SED, we perform a large number (N_MC ≃ 300) of fits of the SED with additional random perturbations. Figure 3 demonstrates the fits of the perturbed SEDs of the 4 pixels of the R9 map. Each model corresponds to one particular set of random perturbations. These perturbations take into account the two main sources of errors, as follows.

• 1.

  The pixel noise at each wavelength (Sect. 4) isassumed to be a normal random independent variable. The noiseis independent from one pixel to the other.

• 2.

  The calibration error is assumed to be a normal random variable, with standard deviation and correlation between wavelengths as described in Sect. 3.4.1. The calibration error from one pixel to the other is correlated.

From a technical point of view, we generate the complete set of independent random variables necessary for the calibration errors, and keep them for our entire analysis. Indeed, one of the advantages of the Monte-Carlo technique is that it allows us to account for complex correlations between variables. As will be demonstrated in Sect. 4.2, the error on the ratio of two quantities depending on the calibration error is often lower than the errors on each individual quantity. It is due to the fact that the calibration error cancels when considering a relative quantity. It is even possible to take this effect into account when comparing the results of two different models, as will be shown in Sect. 4.1.

Figure 4 shows the statistical distributions of the main model parameters, corresponding to the fits of Fig. 3. The parameter value of each Monte-Carlo iteration is the sum of the parameter value of each pixel. The first panel of this figure demonstrates, in particular, that the dust mass has a  ≃50% uncertainty, even for very high signal-to-noise ratio SEDs. In addition it clearly shows the asymmetry of the distribution, which is a result of the non-linearity of the model. The upper end of the statistical distribution of dust masses is less constrained than the lower end, because of the non-linear dependence of the dust mass with the temperature.

Even with a fast running model, computing N_MC ≃ 300 fits for each pixel is CPU intensive, and not necessary. Instead, we use an approximation to reconstruct the probability distributions of each parameter. This method is detailed in Appendix B. It consists of interpolating the pre-computed errors of 30 classes of SEDs, parametrized by their specific power (L_IR/M_dust) and the rms level. Figure 4 compares this reconstruction method to the rigorous Monte-Carlo results. It succeeds in reproducing the correct central value and errors. In particular, it reproduces accurately the skewness of the probability distribution.

3.4.3. Error display

Throughout this paper, each time a numerical quantity is reported as X ≃ a ± b, the two quantities a and b will refer to the mean and standard deviation, or in other words: $\mathit{X}\mathrm{\simeq }\mathrm{⟨}\mathit{X}\mathrm{⟩}\mathrm{\pm }{\mathit{\sigma }}_{\mathit{X}}\mathit{.}$(16)On the other hand, for numerical quantities having a strongly asymmetric distribution, we will quote the error as $\mathit{X}\mathrm{\simeq }{\mathit{a}}_{\mathrm{-}\mathit{b}}^{\mathrm{+}\mathit{c}}$. Only in this case will the quantities refer to the three quartiles, Q₁(X), Q₂(X), Q₃(X), corresponding to values of the repartition function of 1/4, 1/2 and 3/4, respectively: $\mathit{X}\mathrm{\simeq }{\mathit{Q}}_{\mathrm{2}}\mathrm{(}\mathit{X}{\mathrm{)}}_{\mathrm{-}\mathrm{\left[}{\mathit{Q}}_{\mathrm{2}}\mathrm{\right(}\mathit{X}\mathrm{)}\mathrm{-}{\mathit{Q}}_{\mathrm{1}}\mathrm{(}\mathit{X}\mathrm{\left)}\mathrm{\right]}}^{\mathrm{+}\mathrm{\left[}{\mathit{Q}}_{\mathrm{3}}\mathrm{\right(}\mathit{X}\mathrm{)}\mathrm{-}{\mathit{Q}}_{\mathrm{2}}\mathrm{(}\mathit{X}\mathrm{\left)}\mathrm{\right]}}\mathit{.}$(17)The latter error (Eq. (17)) corresponds to a confidence level of 50%, by definition. It is also interesting to consider a higher confidence interval. Most of our error estimates are based on N_MC = 300 Monte-Carlo iterations. It would therefore be meaningless to go down to less than 1/$\sqrt{{\mathit{N}}_{\mathrm{MC}}}\mathrm{\simeq }\mathrm{6}$% error tolerance. We therefore choose to quote the 90% confidence level, defined by the range of the parameter values between 0.05 and 0.95 of the repartition function. From a technical point of view, with N_MC = 300, the limits of this interval are simply the 15th and 286th ordered Monte-Carlo parameter values. We will note this interval: $\mathit{X}\mathrm{\simeq }\mathrm{[}{\mathit{X}}_{\inf }\mathit{,}{\mathit{X}}_{\sup }{\mathrm{]}}_{\mathrm{90}\mathrm{\%}}\mathit{.}$(18)On figures, the 50% error bars will be displayed with a solid line, and the 90% interval with a dashed line.

We note that taking the median (Eq. (17)) gives a central value very close to the best fit, while taking the mean (Eq. (16)) gives a central value systematically shifted from the best fit. Figure 4 demonstrates this effect: the best fit value of the dust mass is 8.62 × 10⁵ M_⊙, very close to the median of the Monte-Carlo iterations (8.76 × 10⁵ M_⊙). On the contrary, the mean, 1.03 × 10⁶ M_⊙, is significantly higher. This is due to the skewness of the probability distribution of the parameter values.

4. The dust mass estimate and its uncertainties

4.1. Systematic discrepancies between the two models

We have applied the model presented in Sect. 3, with the two compositions (“standard” and “AC”), to the maps listed in Table 1. The dust mass spatial distribution obtained with the two models are very similar. Figure 5 shows the distribution of the pixel-to-pixel ratio of the dust masses (“standard” over “AC” models), and starlight intensities, at spatial resolution R4. These distributions are very tight. It shows, that the dust masses obtained with the “AC model” are systematically lower by a factor of $\mathrm{0.3}{\mathrm{8}}_{-0.02}^{\mathrm{+}\mathrm{0.04}}$ than dust masses obtained with the “standard model”, while the starlight intensities are systematically higher by a factor of $\mathrm{2.}{\mathrm{1}}_{-0.2}^{\mathrm{+}\mathrm{0.1}}$. These parameters are tied together. Grains with the “AC model” absorb more light than grains of the “standard model” (Appendix A). They therefore require less mass to account for the observed IR luminosity.

Thus, the results derived with the two compositions give similar trends, but the absolute value of their parameters systematically differ.

4.2. Bias originating in the lack of spatial resolution

In order to quantify the effect of the non-linearity of our SED model, we compute the total dust mass for each map listed in Table 1, with our two grain compositions. Figure 6 shows the resulting trends of dust mass with spatial resolution. These trends are shown for both models. The mass at each resolution is the sum of the masses of all defined pixels. The top panels show the trend of the absolute value of the dust mass. The trends look similar for both models. They appear systematically shifted as discussed in Sect. 4.1. The error bar on each value is large and covers roughly the range of the trend. However, the relative variation of the dust mass (bottom panels of Fig. 6) has significantly smaller error bars. These trends are obtained by normalizing each one of the N_MC Monte-Carlo results of a given spatial resolution, by the corresponding mass for the integrated strip (R10). In that way, most of the calibration error cancels. The only remaining source of uncertainty is the intercalibration error between the various instruments and the rms noise. The dynamics of the trend is unchanged.

The relative trends of Fig. 6 show that measuring the dust mass of an integrated galaxy can give significantly lower values than performing fits of its spatially resolved regions, providing that the spatial resolution is fine enough. For the LMC, the spatially resolved dust mass estimate is  ≃50% higher than the integrated SED fit. This variation is not due to the noise, since the error bars are much smaller than the spread of the trend. Thus, the true mass is probably closer to the high spatial resolution value (R1) than to the integrated flux value (R10). It appears that there is a spatial scale where the trend stabilizes. This transition is an optimal resolution for our model, as it probably provides a correct dust mass, and the noise is lower than at R1. This optimal spatial scale corresponds to R3–R4 (≃27–54 pc; Table 1).

The origin of this trend likely lies in the morphology of the ISM. It could be the result of the dilution of cold massive regions in hotter regions. It is possible that there is a typical spatial scale below which most of the cold regions dominate the SED of the pixels where they lie. It may correspond to the typical scale of molecular complexes. With only a few far-IR/submm constraints, it is difficult to account for these regions when modelling the integrated SED.

To demonstrate the origin of this trend, Fig. 7 compares the pixel distribution of the dust mass surface density, at two spatial resolutions. Taking into account the statistical fluctuations, due to the low pixel number of R7, the two distributions (R1 and R7) agree at low dust mass surface densities (${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{Std}}\lesssim \mathrm{0.06}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}$). However, at high dust mass surface densities (${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{Std}}\gtrsim \mathrm{0.2}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}$), there is a higher fraction of pixels for the high resolution map. This fraction is compensated by an excess of intermediate surface density pixels at low spatial resolution ($\mathrm{0.06}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{pc}}^{-2}\lesssim {\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{Std}}\lesssim \mathrm{0.2}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}$). Thus, the origin of the dust mass underestimation is due to the inability of our model to probe dense regions, at low spatial resolutions. At these spatial resolutions, the mass and average temperature of cold regions are biased by the contribution from hot regions present in the same beam, as the latter are more emissive. On the other hand, at high spatial resolution, cold regions tend to be better separated from hot regions and can therefore be modelled more accurately.

4.3. The gas-to-dust mass ratio crisis: several competitive scenarios

4.3.1. Metal abundance constraints

To be rigorous, we first need to estimate the uncertainty on the metallicity of the LMC. Pagel (2003) compiles the literature for numerous elemental abundances in different regions of the LMC and of the Galaxy. We estimate the error on each element of Table 1 of Pagel (2003), by taking the dispersion of the different measures, excluding cepheids which are very dispersed and not available for all elements. These uncertainties are summarized in Table 5.

We assume that we can reliably derive the metallicity by scaling the mass of O, C and N. In particular, this is supported by the fact that O and C are the major dust constituents. The metallicity of the LMC we adopt is thus: $\begin{array}{ccc}\frac{{\mathit{Z}}_{\mathrm{LMC}}}{{\mathit{Z}}_{\mathrm{\odot }}}& \mathrm{=}& \frac{\mathrm{16}{\left(\frac{\mathrm{O}}{\mathrm{H}}\right)}_{\mathrm{LMC}}\mathrm{+}\mathrm{12}{\left(\frac{\mathrm{C}}{\mathrm{H}}\right)}_{\mathrm{LMC}}\mathrm{+}\mathrm{14}{\left(\frac{\mathrm{N}}{\mathrm{H}}\right)}_{\mathrm{LMC}}}{\mathrm{16}{\left(\frac{\mathrm{O}}{\mathrm{H}}\right)}_{\mathrm{\odot }}\mathrm{+}\mathrm{12}{\left(\frac{\mathrm{C}}{\mathrm{H}}\right)}_{\mathrm{\odot }}\mathrm{+}\mathrm{14}{\left(\frac{\mathrm{N}}{\mathrm{H}}\right)}_{\mathrm{\odot }}}\\ & \mathrm{=}& \mathrm{0.4}{\mathrm{7}}_{-0.07}^{\mathrm{+}\mathrm{0.06}}\mathit{.}\end{array}$(19)This estimate is more accurate than simply scaling the oxygen abundance, as it is usually done.

Let’s check the physical consistency of our dust masses. The Galactic gas-to-dust mass ratio is ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{\odot }}\mathrm{\simeq }\mathrm{158}$ (Zubko et al. 2004). We emphasize that this value is consistent with the elemental depletion patterns. Assuming that the dust-to-gas mass ratio scales with metal abundance (i.e. the dust-to-metal mass ratio is constant), the expected gas-to-dust mass ratio for the LMC is: ${\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathrm{\simeq }\frac{{\mathit{G}}_{\mathrm{dust}}^{\mathrm{\odot }}}{{\mathit{Z}}_{\mathrm{LMC}}\mathit{/}{\mathit{Z}}_{\mathrm{\odot }}}\mathrm{\simeq }\mathrm{33}{\mathrm{9}}_{-41}^{\mathrm{+}\mathrm{55}}\mathit{.}$(20)Assuming that the mass fraction of gaseous heavy elements in the Galaxy is Z_⊙ ≃ 0.017 (Grevesse & Sauval 1998), the solar metallicity dust-to-metal mass ratio is: ${\mathrm{𝒟}}_{\mathrm{\odot }}\mathrm{\simeq }\frac{\mathrm{1}}{{\mathit{G}}_{\mathrm{dust}}^{\mathrm{\odot }}{\mathit{Z}}_{\mathrm{\odot }}}\mathrm{\simeq }\mathrm{0.37.}$(21)It is therefore difficult to understand how the gas-to-dust mass ratio in the LMC could be lower than: ${\mathit{G}}_{\mathrm{dust}}^{\lim \mathit{.}}\mathrm{\simeq }{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathrm{\times }{\mathrm{𝒟}}_{\mathrm{\odot }}\mathrm{\simeq }\mathrm{12}{\mathrm{5}}_{-15}^{\mathrm{+}\mathrm{20}}\mathit{,}$(22)without requiring a larger amount of metals locked-up in grains than what is available in the ISM. These values are summarized in Table 6.

4.3.2. Preliminary: global analysis

Table 7 shows the gas-to-dust mass ratio at each spatial resolution for the two models. These ratios are displayed in Fig. 8. The total gas-to-dust mass ratios given by the “standard model” are too low by a factor of ${\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathit{/}{\mathit{G}}_{\mathrm{dust}}^{\mathrm{Std}}\mathrm{\left(}\mathrm{R}\mathrm{1}\mathrm{\right)}\mathrm{\simeq }\mathrm{3.}{\mathrm{8}}_{-1.0}^{\mathrm{+}\mathrm{1.7}}\mathrm{\simeq }\mathrm{[}\mathrm{2.1}\mathit{,}\mathrm{9.7}{\mathrm{]}}_{\mathrm{90}\mathrm{\%}}$. They are even lower than the physical limit by a factor of ${\mathit{G}}_{\mathrm{dust}}^{\lim \mathit{.}}\mathit{/}{\mathit{G}}_{\mathrm{dust}}^{\mathrm{Std}}\mathrm{\left(}\mathrm{R}\mathrm{1}\mathrm{\right)}\mathrm{\simeq }\mathrm{1.}{\mathrm{4}}_{-0.4}^{\mathrm{+}\mathrm{0.6}}\mathrm{\simeq }\mathrm{[}\mathrm{0.8}\mathit{,}\mathrm{3.6}{\mathrm{]}}_{\mathrm{90}\mathrm{\%}}$. We emphasize here that there is no spatial correlation between the foreground Galactic H i column density (Sect. 2.1) and the gas-to-dust mass ratio deficit of the “standard model”. Therefore, the residual cirrus emission is not responsible for this deficit. Statistically, the “standard model” violates the elemental abundances with a probability of 80%, while the “AC model” is consistent. Therefore, if we assume that our gas mass is correct, then we can conclude that the properties of the “standard model” do not apply to the LMC. This is one scenario.

However, there is a second scenario: the discrepant gas-to-dust mass ratio obtained with the “standard model” could result from the underestimate of the total gas mass. In particular, the mass of molecular gas could have been underestimated. It is known that in low-metallicity environments, the H₂ gas is not properly traced by ¹²CO(J = 1 → 0)_2.6   mm. In these environments, the CO cores are thought to be much smaller relative to their H₂ envelope (H₂ being more efficiently self-shielded). This scenario is supported by the exceptionally high observed [C   ii]_158 μm/¹²CO(J = 1 → 0)_2.6   mm luminosity ratio in dwarf galaxies (e.g. [C   ii]_158 μm/¹²CO(J = 1 → 0)_2.6   mm  ≃ 20   000 in the LMC, compared to [C   ii]_158 μm/¹²CO(J = 1 → 0)_2.6   mm  ≃ 4000 in normal metallicity galaxies; Poglitsch et al. 1995; Israel 1997; Madden et al. 1997; Madden 2000). In principle, the underestimation of the molecular gas mass, using the CO line, could be a factor of  ≃10−100, in these environments (Madden et al. 2011).

Let’s assume that our dust properties are correct, and that ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{Std}}$ is very low due to having underestimated the molecular gas mass. Noting the correction factor accounting for the hypothetical molecular gas not traced by CO and for the fact that our X_CO conversion factor (Sect. 2.3) might be wrong, the total gas-to-dust mass ratio would have to be: $\begin{array}{ccc}\frac{{\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathrm{𝒞}}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}}{{\mathit{M}}_{\mathrm{dust}}^{\mathrm{Std}}\mathrm{\left(}\mathrm{R}\mathrm{1}\mathrm{\right)}}& \mathrm{=}& {\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\\ \mathrm{\Rightarrow }{\mathrm{𝒞}}_{{\mathrm{H}}_{\mathrm{2}}}& \mathrm{\simeq }& \mathrm{10.}{\mathrm{1}}_{-3.6}^{\mathrm{+}\mathrm{5.9}}\mathrm{=}\mathrm{[}\mathrm{4.0}\mathit{,}\mathrm{30.6}{\mathrm{]}}_{\mathrm{90}\mathrm{\%}}\mathit{.}\end{array}$(23)In other words, to explain the discrepant gas-to-dust mass ratio of the “standard model”, we would have to conclude that the molecular gas mass would have been globally underestimated by at least one order of magnitude.

These two alternative scenarios are degenerate, when considering only global values. We therefore need to take into account the redundancy provided by the spatial distribution of the gas-to-dust mass ratio, in order to sort these scenarios out.

4.3.3. Spatial distribution of the gas-to-dust mass ratio

Figure 9 shows the map of IR luminosity for the “standard model”. Since this quantity is the integration of the interpolated observed SED, it depends very little on the model; L_IR distribution of the “standard model” looks almost identical to the “AC model”. That is the reason why we displayed it only for the “standard model”. The upper panels of Fig. 10 show the spatial distribution of the gas-to-dust mass ratio with the two models. A few regions at the two ends of the strip exhibit noisy pixels. They correspond to very low surface densities. In general, there is no particular spatial correlation between the gas-to-dust mass ratio and the CO concentrations. The lower panels of Fig. 10 show the corresponding map of the mass averaged starlight intensity  ⟨ U ⟩ , for both models.

Figure 11 shows the pixel-to-pixel distribution of G_dust for the highest spatial resolution where the gas maps are defined (R4). It appears that the distribution for the “standard model” is shifted to lower values, the pixels are systematically too low, compared to the expected value. In addition, we have built a histogram of the pixels below an arbitrary column density (${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{Std}}\mathit{<}\mathrm{0.1}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}$ and ${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{AC}}\mathit{<}\mathrm{0.03}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}$), in the two lower panels of Fig. 11. We have defined the surface densities at a given spatial resolution by dividing the mass in the pixel by the area of this pixel (Table 1): ${\mathrm{\Sigma }}_{\mathit{X}}\mathrm{=}{\mathit{M}}_{\mathit{X}}\mathit{/}{\mathit{l}}_{\mathrm{pix}}^{\mathrm{2}}$. The two lower panels of Fig. 11 demonstrate that most of the pixels exhibiting high gas-to-dust mass ratios are located in regions of low dust column density, independently of the model used.

Regarding this property, we need to study the variations of the gas-to-dust mass ratio as a function of surface density.

4.3.4. The correlation between gas and dust column densities

In this section, we analyze the correlation between gas and dust using a method similar to that of Planck Collaboration et al. (2011b). It is aimed at identifying the presence of “dark gas” (hereafter DG)¹, as a departure from the linear correlation between gas and dust. We have defined bins of gas mass surface density (${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathrm{\Sigma }}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}$), so that the same number of pixels falls within each bin. We have then computed the average value and the scatter of the dust mass of the pixels within each bin. Figure 12 shows this binned trend on top of the pixel density plot. The binned dust mass surface density is fit with the following function: $\begin{array}{ccc}& {\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{fit}}& \\ & & \mathrm{for}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathrm{𝒞}}_{\mathrm{CO}}^{\mathrm{fit}}{\mathrm{\Sigma }}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\le }{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{DG}}\mathrm{and}{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathrm{𝒞}}_{\mathrm{CO}}^{\mathrm{fit}}{\mathrm{\Sigma }}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{\ge }{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{CO}}\mathit{,}\end{array}$(24)where the free parameters are the following.

• ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{fit}}$ is the fit gas-to-dust mass ratio. It assumes that the actual gas-to-dust mass ratio is the same everywhere in the LMC.

• ${\mathrm{𝒞}}_{\mathrm{CO}}^{\mathrm{fit}}$ is the fit correction factor of the assumed X_CO. Since we have assumed X_CO = 7 × 10²⁰H   cm^-2   (K   km   s^-1)^-1 (Sect. 2.3), then the value of the fit conversion factor is ${\mathit{X}}_{\mathrm{CO}}^{\mathrm{fit}}\mathrm{=}{\mathrm{𝒞}}_{\mathrm{CO}}^{\mathrm{fit}}\hspace{0.17em}\mathrm{7}\mathrm{\times }{\mathrm{10}}^{\mathrm{20}}\mathrm{H}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2} \mathrm{(}\mathrm{K}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}{\mathrm{)}}^{-1}$.

• ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{off}}$ accounts for a possible offset in the gas mass surface density compared to the dust mass surface density.

• ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{DG}}$ is the gas mass surface density above which the dark gas contributes.

• ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{CO}}$ is the gas mass surface density above which the molecular gas is reliably traced by ¹²CO(J = 1 → 0)_2.6   mm, and the dark gas therefore does not contribute anymore. Above this value, the molecular phase is not “dark” anymore.

Equation (24) is not fit between ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{DG}}$ and ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{CO}}$, where the dark gas is assumed to contribute. The total dark gas mass is computed as the difference between the binned trend and the fit of Eq. (24). It is noted ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{DG}}$.

Figure 12 shows the trend and the best fit for each model. Although the departure between the line fit and the binned trend, in the $\mathrm{\left[}{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{DG}}\mathit{,}{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{CO}}\mathrm{\right]}$ is visible, its deviation is lower than the typical dispersion of the correlation. However, the parameter values and their uncertainties, given in Table 8, show this departure is significant. It appears that the gas-to-dust mass ratios derived from these simple fits are consistent with the pixel-to-pixel values (Table 7). Therefore, this analysis tends to confirm that the “standard model” violates elemental abundances in the LMC, while the “AC model” is physically valid. The fact that the parameter ${\mathrm{𝒞}}_{\mathrm{CO}}^{\mathrm{fit}}$ is consistent with unity indicates that our adopted X_CO conversion factor (X_CO = 7 × 10²⁰H   cm^-2   (K   km   s^-1)^-1; Sect. 2.3) is probably not incorrect. Finally, the dark gas mass fraction, around 10%, seems to be moderate.

There are several hypotheses entering into Eq. (24). First, its interpretation implicitly relies on the translation between Σ_dust and A_V, the extinction magnitude in V band, because it defines the limits on the dark gas regime (${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{DG}}$ and ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{CO}}$; Eq. (24)). However, at the spatial scales considered here (54 pc), molecular clouds are not resolved, and several phases are mixed within each pixel. Therefore, the dust mass surface density is in principle a biased estimator of the A_V. The value of A_V derived from the dust mass surface density, assuming a uniform dust distribution, is actually always going to be lower than the A_V of a molecular cloud that would lie in the pixel. Second, the approach of Eq. (24) is a perturbative method. It is correct only if the fraction of dark gas is small. In particular, this method would fail if dark gas was present at low Σ_gas. Finally, it relies on the assumption of a uniform gas-to-dust mass ratio throughout the entire galaxy. In fact, our G_dust spatial distributions (Fig. 10) are not noise maps, they contain clear structures. Moreover, our error analysis demonstrates that the amplitude of these structures is larger than the typical error bar on an individual pixel value. And the possible offset in background subtraction between the gas and dust maps (estimated in Table 8) is small enough to affect only the low surface brightness pixels. Consequently, this prompts us to further scrutinize the observed variations of G_dust as a function of the physical conditions.

4.3.5. Variations of the gas-to-dust mass ratio with physical conditions

In this section, we describe the observed trends of the gas-to-dust mass ratio as a function of several parameters. Figure 13 shows the pixel-to-pixel variations of the observed gas-to-dust mass ratio, as a function of several tracers of the physical conditions, for the two models. We note that the trends are similar with the two models, but the quantities proportional to the dust mass appear to be scaled down by a factor of  ≃0.38 with the “AC model”. On the contrary,  ⟨ U ⟩  appears to be scaled up by a factor of  ≃2.1 (see also Fig. 5).

The two left panels of Fig. 13 show the variation of G_dust as a function of the gas mass surface density Σ_gas. There is a large scatter within each bin of Σ_gas, and the scatter is higher at low surface densities, since the signal-to-noise ratio is lower and offsets in the zero level of the gas and dust masses can impact the ratio. The gas-to-dust mass ratio at high gas mass column densities (Σ_gas ≳ 10 M_⊙   pc^-2) appears to deviate slighty from a constant value, as discussed in Sect. 4.3.4. However, this deviation is smaller than the dispersion. Overall, this trend does not show any general correlation between the gas-to-dust mass ratio and the gas mass surface density. The central values of G_dust are on average below the limit value with the “standard model”, but are found around the expected value with the “AC model”.

The two central panels of Fig. 13 show the variation of G_dust with the dust mass surface density Σ_dust. There is a clear anticorrelation between the two quantities. Although, the trend at low surface densities (Σ_dust ≲ 0.01 M_⊙   pc^-2) is very scattered and sensitive to the offset ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{off}}$ (Table 8), the rest of the trend is clear. The fact that the gas-to-dust mass ratio correlates better with the dust mass column density than with the gas mass column density is the sign that the origin of the variations of G_dust is not directly linked with the observed gas content. Most of the trend is below the limit value, with the “standard model”. On the contrary the values for the “AC model” are all realistic.

Finally, the two right panels of Fig. 13 show the variation of G_dust with the mass averaged starlight intensity  ⟨ U ⟩ . There is a clear relation between the two quantities. Contrary to the central panels (G_dust vs. Σ_dust), these two quantities are independent. Moreover,  ⟨ U ⟩  is a specific quantity, while Σ_dust is integrated over the line of sight. It appears that  ⟨ U ⟩  is the parameter giving the best trend with G_dust, exhibiting larger dynamics and smaller scatter than the trends with Σ_gas and Σ_dust. We note that, with both models, the slope of the correlation (in bi-logarithmic representation) is not constant. Schematically, for  ⟨ U ⟩ ^Std ≲ 1 (or  ⟨ U ⟩ ^AC ≲ 2), G_dust ∝  ⟨ U ⟩ . Then, for 1 ≲  ⟨ U ⟩ ^Std ≲ 4 (or 2 ≲  ⟨ U ⟩ ^AC ≲ 8), the trend flattens. There is no significant variation in this range. At larger values, G_dust ∝  ⟨ U ⟩  again, but the pixels are very dispersed in this range, and the statistics are limited.

In summary, the observed variations of G_dust do not appear to be significantly traced by the observed gas content. On the contrary, these variations are well correlated with tracers of the dust physical conditions, in particular their irradiation conditions. It suggests that the origin of these variations is linked in some way with the radiation energy density. We are now going to explore the possible interpretations of these variations.

4.3.6. Physical interpretation: dark gas, modified submm opacities or enhanced dust condensation?

To interpret the trends of Fig. 13, we need to list the possible physical processes responsible for variations of the observed gas-to-dust mass ratio, as well as the known biases affecting these trends. The known biases of our trends are the following.

• Bias 1.

  There is a potential offset between the gas and dust masszero level, due to differences in background subtractions. This offset has an impact on low surface density pixels (or high ⟨ U ⟩ ). It is quantified in Table 8.

• Bias 2.

  There is a known bias of our SED fitting. At high average starlight intensities ( ⟨ U ⟩  ≳ 10), our model tends to underestimate the dust mass. This effect is demonstrated in Fig. B.2 of Appendix B. An artificial correlation between  ⟨ U ⟩  and G_dust is expected in this regime.

The possible causes of variations of the observed gas-to-dust mass ratio are the following.

• Dark gas.

  As discussed in Sect. 4.3.2, lowG_dust could be the sign of massive “dark gas” reservoir. As shown in Sect. 4.3.5, the fact that clear variations of G_dust are correlated with Σ_dust, but not with Σ_gas, suggests that the H i and CO traced gas is blind to actual variations of G_dust. Assuming that our X_CO is correct, the decrease of G_dust with Σ_dust could be partly due to the presence of a dark gas component more massive than what was determined in Sect. 4.3.4. The latter could have been underestimated because it was lost in the scatter of the relation.

• Dust-to-metal ratio.

  Variations of G_dust could be due to local variations of the dust-to-metal mass ratio (Eq.  (21 )). In principle, the gas-to-dust mass ratio could go down to ${\mathit{G}}_{\mathrm{dust}}^{\lim \mathit{.}}$ (Eq. (22)) and go up to infinity. Assuming an efficient mixing of freshly formed dust grains throughout the ISM, an increase of G_dust would be the sign of an enhanced destruction, mainly by SN blastwaves. This is probably the case at low surface density, although it is degenerate with biases 1 and 2. On the other hand, a decrease of G_dust would then be the sign of an increase of the condensation efficiency of metals onto dust. This effect is degenerate with the dark gas component. Assuming there is no dark gas, for this scenario to be realistic, the metal condensation onto dust is expected to be enhanced in regions where the gas density is higher. However, the left panels of Fig.  13 show that there is no clear correlation between G _dust and the gas mass surface denstity. This effect should therefore be minor. Although the dust-to-metal mass ratio may vary in the LMC, our data indicate that this potential variation is not at the origin of the trend, except maybe for low surface density.

• Emissivity.

  The variations of G_dust could be due to variations of the emissivity of the grains in different regions. This emissivity variation would lead us to systematically misestimate the dust mass in regions where it would differ from the assumed opacity. These variations are poorly known and more difficult to constrain. The most well-known of these processes, grain-grain coagulation, leads to an increase of the far-IR opacity (e.g. Stepnik et al. 2003). It happens preferentially in dense clouds. However, the left panels of Fig. 13 do not show a clear correlation between G_dust and Σ_gas. Therefore, this effect does not seem to be at the origin of the overall variations of G_dust. Another effect is the increase of the emissivity with grain temperature (e.g. Meny et al. 2007). This effect would lead us to overestimate the dust mass in regions where the dust is hot. It would therefore lead to an anticorrelation between the observed G_dust and  ⟨ U ⟩ . Therefore, this effect can not be at the origin of our trend. Although this review is not complete, we will assume that local variations of the grain emissivity are not responsible for the observed variations of G_dust.

Considering these previous effects, we now make several simple assumptions to interpret the trends of G_dust with  ⟨ U ⟩  (right panels of Fig. 13). We assume that we can classify our pixels within the three following regimes, characterized by their average starlight intensity.

• A.

  We assume that there are embedded regions, characterized bylow starlight intensities(1 ≲  ⟨ U ⟩ ), where the observed gas-to-dust mass ratio is underestimated due to the presence of an undetected gas phase. The proportion of this overlooked gas rises, when  ⟨ U ⟩  decreases.

• B.

  We assume that there are regions, characterized by moderate starlight intensities (1 ≲  ⟨ U ⟩  ≲ 10), corresponding to the H i dominated regime, where there is no dark gas. The observed gas-to-dust mass ratio is assumed to be correct in these regions, and can be used as a reference.

• C.

  We assume that the gas-to-dust mass ratio may be higher in very diffuse regions, where shocks might have recently destroyed the grains. In these regions (⟨ U ⟩  ≳ 10), the biases of the modelling can also be significant, as demonstrated in Appendix B and Fig. B.2. An artificial correlation between  ⟨ U ⟩  and G_dust is expected.

Figure 14 illustrates this scenario. The two panels show the binned trends of Fig. 13 with the starlight intensity, for the two models. The grey circles with error bars show the original trend, while the black circles with error bars show the trend corrected for the offset ${\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathrm{off}}$ derived in Table 8. We fit the latter. We decompose the trends into the three zones described above, by fitting the following function: $\log {\mathit{G}}_{\mathrm{dust}}^{\mathrm{fit}}\mathrm{=}\{\begin{array}{c}\\ {\mathit{a}}_{\mathrm{dark}}\log \frac{\mathrm{⟨}\mathit{U}\mathrm{⟩}}{\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{dark}}}\mathrm{+}\log {\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}\mathrm{for}\frac{\mathrm{⟨}\mathit{U}\mathrm{⟩}}{\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{dark}}}\mathit{<}\mathrm{1}\\ \\ \log {\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}\mathrm{for}\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{dark}}\mathrm{\le }\mathrm{⟨}\mathit{U}\mathrm{⟩}\mathrm{\le }\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{diff}}\\ \\ {\mathit{a}}_{\mathrm{diff}}\log \frac{\mathrm{⟨}\mathit{U}\mathrm{⟩}}{\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{diff}}}\mathrm{+}\log {\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}\mathrm{for}\frac{\mathrm{⟨}\mathit{U}\mathrm{⟩}}{\mathrm{⟨}\mathit{U}{\mathrm{⟩}}_{\mathrm{diff}}}\mathit{>}\mathrm{1}\mathit{,}\end{array}$(25)where a_dark, a_diff, ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}$,  ⟨ U ⟩ _dark and  ⟨ U ⟩ _diff are free parameters. ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}$ is the “reference” gas-to-dust mass ratio, i.e. the gas-to-dust mass ratio in the diffuse ISM, and  ⟨ U ⟩ _dark and  ⟨ U ⟩ _diff are the starlight intensities defining these zones (Fig. 14). Table 9 gives the main parameters of these fits for each model.

This point of view leads to the same conclusion than in Sect. 4.3.4, namely that the “standard model” is unphysical. Indeed, the derived reference gas-to-dust mass ratio, ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}\mathrm{\simeq }\mathrm{11}{\mathrm{4}}_{-14}^{\mathrm{+}\mathrm{25}}$ (Table 9), in regions where dark gas is unlikely, is on average below the hard limit, ${\mathit{G}}_{\mathrm{dust}}^{\lim \mathit{.}}\mathrm{\simeq }\mathrm{125}$. On the contrary, the “AC model” gives a reference ${\mathit{G}}_{\mathrm{dust}}^{\mathrm{ref}}\mathrm{\simeq }\mathrm{30}{\mathrm{6}}_{-38}^{\mathrm{+}\mathrm{34}}$ close to the expected value, based on the metallicity of the LMC, ${\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathrm{\simeq }\mathrm{339}$. More importantly, this reference ratio does not violate the elemental abundances, with the “AC model”.

In summary, now that we have taken into account the various competing processes, we can safely conclude that the “standard model”, which works for the Milky Way, does not apply to the LMC. It means that the chemical composition of the grains in the LMC is systematically different than that of the Galaxy. The LMC grains have on average a larger submm opacity. Although the possible compositions are numerous, this work gives a plausible one, respecting the elemental abundances, based on optical properties of grains that have been observed in the laboratory (ACAR amorphous carbons in lieu of graphite; Zubko et al. 1996).

4.3.7. Remarks on the submillimiter emissivity index

First, as detailed in Appendix A, our “AC model” has a lower far-IR emissivity index (β ≃ 1.7) than the “standard model” (β = 2). In our case, β does not vary much with wavelength in the Herschel range. However, in general, β is a function of the wavelength, and there are grain species showing strong variations in this range (e.g. Quinten et al. 2002).

Second, the emissivity of our grain model is what we call the intrinsic emissivity. On the other hand, the isothermal fit of an SED with free emissivity index provides the effective emissivity index (noted β_eff. by us). The two are similar (β = β_eff.) only if the isothermal approximation is valid. On the contrary, if the spatial resolution is such that there is likely a mix of temperatures within the beam, then the two values are going to differ. The SED will be broadened by the temperature distribution and β_eff. will likely be lower than β. In the LMC, we are mostly in the second case. Due to this degeneracy, the shapes of the SED do not provide direct constraints on β. Other constraints have to be invoked. In our case, we used the constraints on the gas-to-dust mass ratio to discriminate between different β.

We note that Aguirre et al. (2003), using COBE observations of the LMC, reached the same qualitative conclusion as derived here. They concluded that the grain opacities were different in the LMC than in the Milky Way, although they adopted a different set of grain species. However, these authors based their analysis only on isothermal fits, and did not use the gas-to-dust mass ratio as a constraint. Therefore their demonstration was incomplete. Similarly, Gordon et al. (2010, using Herschel and Planck Collaboration et al. (2011a, using Planck showed that the effective β of the LMC was around 1.5. However, they did not constrain the intrinsic β.

Thus, in light of these precisions, we have shown, for the first time, that the intrinsic far-IR/submm emissivity index of the grains in the LMC is lower than in the Galaxy.

4.3.8. Constraints on the dark gas content

Using the point of view of Fig. 14, we can put an upper limit on the dark gas content. Since we have shown in Sect. 4.3.6 that most of the variation of the gas-to-dust mass ratio could originate in an overlooked gas component, we can estimate the mass of this component, from this point of view.

To provide an optimistic estimate, we define the mass of overlooked gas ${\mathit{M}}_{\mathrm{gas}}^{\mathit{/}\mathrm{DGup}}$ by the mass added to the observed gas mass ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}$, in order to match the expected gas-to-dust mass ratio ${\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}$ in regions where the uncorrected G_dust is lower than ${\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}$: ${\mathit{M}}_{\mathrm{gas}}^{\mathrm{DGup}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{(}{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathrm{-}{\mathit{G}}_{\mathrm{dust}}\mathrm{)}\mathrm{\times }{\mathit{M}}_{\mathrm{dust}}& \mathrm{if}{\mathit{G}}_{\mathrm{dust}}\mathrm{\le }{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\\ & \\ \mathrm{0}& \mathrm{if}{\mathit{G}}_{\mathrm{dust}}\mathit{>}{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}\mathit{.}\end{array}$(26)Technically, this dark gas includes the CO-free H₂, as well as a fraction of the CO-associated H₂ where the X_CO factor is higher than our adopted value. This definition implies that this dark component covers most of the strip with the “standard model”, while it is going to cover only a fraction of it with the “AC model”. Table 10 gives the mass fraction of this component, for both models, at each spatial resolution. The most reliable value is probably the value at R4 of the “AC model”.

Therefore, comparing the conservative estimate of Table 8 and the optimistic estimate of Table 10, with their respective uncertainties, we can bracket the dark gas mass fraction ${\mathit{f}}_{\mathrm{gas}}^{\mathrm{DG}}\mathrm{=}{\mathit{M}}_{\mathrm{gas}}^{\mathrm{DG}}\mathit{/}\mathrm{(}{\mathit{M}}_{\mathrm{gas}}^{\mathrm{H}\hspace{0.17em}\mathrm{i}}\mathrm{+}{\mathit{M}}_{\mathrm{gas}}^{{\mathrm{H}}_{\mathrm{2}}}\mathrm{)}\mathrm{\simeq }\mathrm{10}\mathrm{-}\mathrm{100}\mathrm{\%}$ in the strip, with the “AC model”.

For comparison, Bernard et al. (2008) found that this fraction was 182%. The reason for this discrepancy might be due to the fact that the region we have modelled here is relatively poor in dark gas, compared to the entire LMC.

4.4. What the SPIRE_500 μm excess is not

As mentioned in Sect. 3.3, we have not fitted the SPIRE_500 μm excess present in our data. We have excluded this waveband, in order to avoid being biased by this effect, since its origin is still unknown. This excess can, a priori, affect our previous result concerning the grain properties, for two different reasons. First, this excess could originate in colder dust that our model has not accounted for. Our dust mass would then have been biased. Second, if the origin of this excess is an unkown mechanism, not related to cold dust, it could also affect shorter wavelength bands, and bias our dust mass.

Although SPIRE_500 μm was not used as a constraint, for our dust SED model, it is possible to study the behaviour of its excess with the physical conditions.

4.4.1. Systematic analysis of the excess

Let’s define the absolute and relative excesses of the SPIRE_500 μm band, respectively by: $\begin{array}{ccc}{\mathrm{\Delta }}_{\mathrm{500}}& \mathrm{=}& \\ {\mathit{r}}_{\mathrm{500}}& \mathrm{=}& \end{array}$Figure 15 shows four SEDs with different values of the excess.

Table 11 shows the value of the excesses, as a function of the spatial resolution. In general, the excess tends to be slightly smaller with the “AC model”, since it has a slightly flatter submm slope. The relative excess tends to rise slightly when the spatial resolution decreases, but this trend is not statistically significant. We emphasize that r₅₀₀ being a relative quantity, most of the calibration error cancels, and the errors quoted in Table 11 are uncorrelated between different spatial resolution. Most of the error in this excess comes from the noise in the diffuse emission. The most reliable estimates of the excess are therefore those of R10 (integrated strip): ${\mathit{r}}_{\mathrm{500}}^{\mathrm{Std}}\mathrm{=}\mathrm{16.}{\mathrm{1}}_{-7.8}^{\mathrm{+}\mathrm{6.4}}\mathrm{\%}$ and ${\mathit{r}}_{\mathrm{500}}^{\mathrm{AC}}\mathrm{=}\mathrm{14.}{\mathrm{1}}_{-6.9}^{\mathrm{+}\mathrm{6.1}}\mathrm{\%}$. This excess is only marginally detected in the global SED. Since, the excess is relatively similar for both models, and that we have previously shown that the “AC model” is more realistic, we will discuss the excess only for the “AC model”, in the rest of this section.

Figure 16 compares the spatial distribution of the absolute SPIRE_500 μm excess (Eq. (27), in colors), with the relative excess (Eq. (28), with contours). We first note that the absolute excess is not homogeneous, and contains structures associated with the known features of the LMC, in particular N   44, and the north superbubble. We also note that the absolute value of the excess is much larger than the typical noise (Table 4). We emphasize that the map has been background subtracted in such a way that the fluxes on the north and south of the strip are zero. These arguments show that this excess can not be attributed to emission external to the LMC, like CMB fluctuations or cirrus foreground. The distribution of the relative excess r₅₀₀ also shows structures, that appear to be roughly anticorrelated with the IR luminosity.

Let’s look in more details at the correlation of the excess with different quantities, in order to decipher its origin. Figure 17 shows the relation of the absolute excess to the gas phases. It shows that the absolute excess is not correlated with the molecular gas, but is correlated with the H i column density. It therefore suggests that the SPIRE_500 μm excess is associated with the atomic medium but not with dense phases. To go further, we need to look into the various interpretations that have been proposed, for submm excesses in galaxies.

4.4.2. Consistency test of different interpretations

A submm emission excess has been reported in different systems. We warn the reader that this excess extends up to cm wavelengths. It is still uncertain if the same physical process is responsible for the entire wavelength range, or if it is the combination of several phenomena.

Reach et al. (1995) first reported a long wavelength excess in the COBE/FIRAS spectrum of the Galaxy. This excess could be fit with a very cold dust component (4–7 K). However, the authors rejected this solution, since it was located at high Galactic latitudes. This excess was fit by Li & Draine (2001) invoking a long wavelength enhancement of the opacity of their silicate grains (see also Appendix A). Galliano et al. (2003) discovered a submm excess in the SCUBA and MAMBO observations of the blue compact galaxy NGC1569. Then several studies reporting similar excesses in galaxies, mainly low-metallicity dwarf galaxies were published, e.g.: Dumke et al. (2004), Galliano et al. (2005), Bendo et al. (2006), Galametz et al. (2009, 2010), O’Halloran et al. (2010). The relative intensity of this excess appears to decrease with metallicity (e.g. the sequence SMC, LMC, MW). Planck Collaboration et al. (2011a) presented the global excess at long wavelength.

The main interpretations of this excess, found in the literature, are the following.

• 1.

  Gallianoet al. (2003, 2005)discussed the consequences of a very cold dust component,showing that it was not inconsistent if the emissivity index of thesegrains was β ≃ 1. This scenario would imply a small number (at most a few hundreds) of very dense, parsec-size clumps, containing a large fraction (40−80%) of the ISM mass.

• 2.

  Meny et al. (2007) showed that the excess was successfully fit with a physical model of temperature dependent grain emissivity. Their model implies that the emissivity index decreases, when the temperature of the grain increases.

• 3.

  Bot et al. (2010) have shown that the long wavelength spectrum of the excess renders the sole very cold hypothesis impossible. It has to be combined with another process to be realistic: e.g. very cold dust and spinning dust grains (Draine & Lazarian 1998). Alternatively, Planck Collaboration et al. (2011a) proposed that the excess in the SMC could be a combination of the Meny et al. (2007) grain model and of spinning grains. However, Planck Collaboration et al. (2011a) attributed the excess in the LMC to CMB fluctuations. Bot et al. (2010) explained the excess in the LMC and SMC with a combination of 12 K dust and spinning grains (Ysard & Verstraete 2010).

Here, we discuss r₅₀₀ in light of these scenarios.

Since the actual value of this excess is very sensitive to the noise, let’s zoom into the brightest star forming region of our observations: N   44 and its surroundings. Figure 18 shows the central distribution of the excess around N   44, compared to the starlight intensity and the dust mass column density. Figure 19 shows the corresponding pixel-to-pixel variations. It also demonstrates the typical error bars on individual pixels, showing that the excess is significant when looking at spatially resolved regions. Both  ⟨ U ⟩  and Σ_dust are correlated with r₅₀₀. However, the correlation is better with Σ_dust. This is confirmed by the spatial distribution: although, on average the excess corresponds to high starlight intensity regions, there is one spot north of N   44, with no excess and high starlight intensities. Therefore, our observations indicate that the excess is primarily associated with diffuse regions. As a consequence, this excess can not come from contamination by the ¹²CO(J = 5 → 4)_520   μm line.

In regard to these trends, the very cold dust hypothesis is very unlikely. Indeed, we would expect very dense clumps to be associated with dense regions. On the contrary, the observations show that most of these clumps would be present in the diffuse ISM, and very few or none of them in dense regions. However, this is only a qualitative argument showing that very cold dust is unlikely. To rule out this component, we develop the following, more quantitative point of view. In order to reach temperatures below 10 K, very cold dust has to be efficiently shielded. And, the dust responsible for the shielding is going to absorb the ambient radiation and reach warm temperatures, emitting in the wavelength range that has been modelled in the previous sections. We can therefore estimate the mass of shielding dust necessary to allow the very cold dust to reach low temperatures. By comparing the order of magnitude of the minimum mass of shielding dust to the observed mass in each pixel, Appendix D shows that very cold dust can reasonably be ruled out.

Incidentally, by showing: (i) that very cold dust is unlikely to contribute significantly to SPIRE_500 μm and (ii) that the SPIRE_500 μm excess is found mainly in diffuse regions, we have demonstrated that this excess does not have any impact on our main discussion about dust mass (Sect. 4.3).

5. Summary and conclusion

In this paper, we have presented the modelling of the spatially resolved Spitzer/IRAC, Spitzer/MIPS and Herschel/SPIRE data of a strip covering one quarter of the LMC. The purposes of this work was to:(i) systematically study all the effects leading to inaccuracies or biases affecting the dust mass estimate of a galaxy, in the Herschel era; (ii) explore the peculiar ISM properties of the LMC. Our main results are the following.

• 1.

  We have presented an empricial model to fit IR/submm SEDs.This model adopts realistic grain properties, and accounts for apossible distribution of starlight intensity in the region studied.We have described in detail the general propagation of theobservational errors (noise and calibration) through the SEDfitting process. We have shown that, even when thesignal-to-noise ratio is high, the errors on the dust mass areimportant and strongly asymmetric (typically${\mathrm{\simeq }}_{\mathrm{-}\mathrm{25}\mathrm{\%}}^{\mathrm{+}\mathrm{40}\mathrm{\%}}$). However, we have shown that relative quantities (ratio of two parameters) can have small error bars (typically ${\mathrm{\simeq }}_{\mathrm{-}\mathrm{7}\mathrm{\%}}^{\mathrm{+}\mathrm{10}\mathrm{\%}}$) due to the partial cancelation of the correlated calibration errors.

• 2.

  By modelling the same set of maps, but with different pixel sizes, we have shown that the lack of spatial resolution can lead to a systematic underestimate of the dust mass by  ≃50%. This bias could be the result of the veiling of cold components by the emission from warmer regions. Although the amplitude of this bias is probably specific to the type of environment found in the LMC and to the wavelength coverage of our data set, we believe that this trend with spatial resolution is general. Modelling the integrated SED of a galaxy leads to underestimating the dust mass.

• 3.

  We have performed our SED fitting with two sets of grain composition, in order to explore the sensitivity of the dust mass estimate to the submm emissivity of the grains. We have shown that both compositions give equally good fits, but that the “standard” grain composition (graphite & silicate), that works for the Milky Way, fails to give physically realistic results in the LMC. It violates the elemental abundances. We discuss the fact that this discrepancy could either be due to modified submillimeter grain opacities (our AC model: amorphous carbon & silicate), or to the presence of an unaccounted for gas reservoir. In particular, we have shown that there is a degeneracy between the intrinsic grain emissivity and their temperature distribution, and that considerations on the mass are a powerful way to remove some of this degeneracy.

• 4.

  The detailed spatial analysis of the variations of the observed gas-to-dust mass ratio has proven the need for grain opacities different than those of the Galaxy. More precisely, our analysis has demonstrated that grains in the LMC have on average a larger far-IR/submm opacity. We propose a physical dust model that is consistent with these properties (emissivity index β ≃ 1.7, and κ_abs(160 μm) = 1.6 m² kg^-1), although we insist that this particular composition is not a unique solution.

• 5.

  We have shown that the mass averaged starlight intensity is a better tracer of regions where the observed gas-to-dust mass ratio is depleted. Comparing different methods, we have constrained the amount of “dark gas” (i.e. unaccounted for by [H i]_21 cm, and ¹²CO(J = 1 → 0)_2.6   mm). This reservoir accounts for  ≃ 10–100% of the gas mass of the strip. Simply correcting the standard CO-to-H₂ conversion factor is not sufficient to account for this dark gas reservoir, since it is abundant in regions where no CO detection is reported.

• 6.

  We have analyzed the excess emission in the SPIRE_500 μm waveband, which has been reported before in the LMC and other low-metallicity galaxies. On average the amplitude of the relative excess is  ≃15%, and it can vary spatially between 0 and  ≃40%. By scrutinizing its behaviour, we show that this excess is mainly anticorrelated with the density of the ISM. We show that this excess can not be due to CMB fluctuations, nor to massive amounts of very cold dust. However, the nature of this excess remains unknown.

Considering the uncertainties in grain optical properties, our study demonstrates that, when modelling the IR/submm SED of a galaxy, the derived dust mass can be considered as a constraint rather than a result. Indeed, we have shown that standard Galactic grain properties were leading to unphysical masses. We therefore have selected another model which was giving reasonable results. However, going from this conclusion to providing the reader with an actual dust mass or gas-to-dust mass ratio for the LMC would be a circular process. It is now clear that grain properties can vary significantly from a galaxy to another. The range of gas-to-dust mass ratios that results from using different grain species makes the estimate of the absolute dust mass irrelevant, without independent constraints. Such independent constraints could be the elemental depletions, or spectroscopic information on the precise chemical composition of the grains, or simultaneous extinction measurements.

Appendix A: The grain properties of our models

A.1. The submillimetre opacities

The top panel of Fig. A.1 compares the absorption opacities of the grain mixtures of our two models. These opacities are the sum of PAHs, carbon grains and silicates. The effective submillimeter emissivity index β is defined by the logarithmic index of the opacity: ${\mathit{\kappa }}_{\mathrm{abs}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{\propto }{\mathit{\lambda }}^{\mathrm{-}\mathit{\beta }}\mathit{.}$(A.1)It is β ≃ 2 for the “standard model” and β ≃ 1.7 for the “AC model”, although β for amorphous carbons alone is even lower. Table A.1 gives the submilimeter properties of our two models approximated by Eq. (A.1).

The bottom panel compares the infrared emission of the compositions, for the same starlight intensity. It demonstrates that the “AC model” has more emissivity, especially in the submm. It therefore allows us to fit the same observed SED with less mass, and slightly hotter grains.

Fitting the COBE/FIRAS high latitude Galactic cirrus, Li & Draine (2001) had to modify the imaginary part of the dielectric function of their silicate grains, at wavelengths greater than 250 μm. This modification would have a very limited effect on our conclusion, since our constraints go up to 350 μm only, where this increase is very limited (it is  ± 12% in the 250 μm ≤ λ ≤ 1000 μm range). This modification actually lowers the emissivity between 250 and 850  μm, and increases it, at λ > 850 μm. For our dust mixture, the amplitude of this modification is even lower, since the contribution of silicates to the far-IR is lower than for the Li & Draine (2001) model. The difference between the “standard model” and the one using the modified silicate of Li & Draine (2001) is invisible in Fig. A.1. The physical origin of this excess may be similar to the one we find here. However, its amplitude is much larger in the LMC, and manifests at shorter wavelengths (λ ≳ 100 μm).

A.2. Size distribution considerations

It is important to note that the size distribution used here does not include grains larger than 0.35 μm, and the contribution to the emission of grains larger than a ≳ 0.1 μm is negligible. Large grains tend to have a flat UV-visible opacity, and therefore, to have a lower equilibrium temperature than smaller equilibrum grains, exposed to the same ISRF. Therefore, adding larger grains would increase the submillimetre emissivity of the model, and allow us to fit the SPIRE fluxes, without having to go to low starlight intensities. However, having large grains or having colder dust would give the same discrepancies in terms of gas-to-dust mass ratios. Moreover, a significant excess of large grains would flatten the UV rise of the extinction curve, contradicting the observations of several lines of sight within the LMC (Gordon et al. 2003).

Another feature of our model is that we have been forced to lower the abundance of non-PAH small grains (a < 10 nm; both carbon and silicate grains) by a factor of 2. Without this modification, we were not able to get good fits of the MIPS_{24  μm} of the diffuse regions. This modification has a very minor effect on the dust mass (≲ 10%), and it is systematic, thus it has no impact on our conclusions. However, this is puzzling since the fit of the extinction curves of the LMC indicates a larger fraction of small grains (e.g. Weingartner & Draine 2001). The fact is that the Zubko et al. (2004, BARE-GR-S) model has a higher small grain contribution in MIPS_24 μm than the Draine & Li (2007). On the other hand, replacing graphite by amorphous carbon, as shown in this paper, allows us to get rid of the 30 μm graphite feature, and decrease the contribution of small grains in the MIPS_24 μm band, making the 24 μm fit better without having to alter the size distribution. This is another indirect consistency check of the conclusion of this paper.

A.3. Overview of the derived PAH properties

Although the scope of our paper was not to discuss the PAH properties, their abundance is a natural output of our model. In this section, we summarize these results.

Table A.2 shows the PAH mass fractions for the two models. This parameter is relatively well constrained and does not vary significantly with spatial resolution. Indeed, it depends essentially on the IRAC_8 μm-to-total-IR luminosity ratio. The mass fraction for the “AC model” is systematically higher, since the bulk of the dust is more emissive.

Figure A.2 shows the spatial distribution of the mass fraction of PAHs, f_PAH (Eq. (3)). We confirm the Paradis et al. (2009) results showing an excess of PAH emission toward the stellar bar. However, f_PAH is biased by the fact that, in absence of detailed mid-IR spectrum, we have arbitrarily fixed the charge fraction to 1/2. The charge fraction controls the emissivity of the C-C modes (Galliano et al. 2008b). Therefore, it is difficult to uniquely interpret this excess emissivity by a local increase of the PAH abundance.

Appendix B: Details concerning the error analysis

B.1. Classes of observed SEDs

The error analysis presented in Sect. 3.4 would imply having to fit the 156   577 pixels (Table 1) 300 times with the two models, leading to a total of  ≃10⁸ fits. However, all these SEDs have a lot of similarities, and it is not necessary to perform the Monte Carlo iterations on each pixels. Here we describe the approximation method we demonstrated in Fig. 4.

We start from the actual fit of the unperturbed SED of each pixel of Table 1. We order these SEDs according to their specific emitted power , making 30 logarithmic bins. Figure B.1 shows the classes of observed SEDs, for the “standard model”. At each waveband λ₀, the central value is the median of the normalized monochromatic power λ₀L_λ0/L_IR of each pixel within the considered bin of specific power. The error shows the dispersion of the pixels within the bin. On first approximation, the specific power is proportional to the mass averaged starlight intensity  ⟨ U ⟩ . That is the reason why the SEDs of Fig. B.1 are nicely ordered according to the wavelength peak of the dust emission. Notice that the dispersion at most wavebands is small compared to the variation spread by the different classes. The only wavelength range where it is not true is the near-IR (λ ≲ 5 μm), where the stellar emission dominates.

B.2. Class interpolation and error estimate

We perform the Monte Carlo error analysis of Sect. 3.4.2 on each SED class, for 6 different noise levels, spanning the whole range of observed signal-to-noise ratios, S/N. We perform this study with the two models. Then, for each pixel of Table 1, the errors on the parameters are determined by logarithmically interpolating the precomputed errors in and S/N.

The validity of this interpolation method is demonstrated in Fig. 4. We show that this method reproduces well the central value and error bars of each parameter, since it reproduces well the skewness of the probability distribution.

B.3. Biases of the SED fits

The large database of Monte-Carlo fits of the SEDs of Fig. B.1 can be used to study the systematic effects of our method. Figure B.2 shows the bias on the dust mass, as a function of the mean equilibrium grain temperature.

The two panels of this figure correspond to the two models. The x-axis is the “actual”  ⟨ U ⟩ , i.e. the value of  ⟨ U ⟩  derived from the unperturbed SEDs of Fig. B.1. Each value of  ⟨ U ⟩  defines one of the 30 classes, since . The y-axes show the ratio between: ${\mathit{M}}_{\mathrm{dust}}^{\mathrm{measured}}$, which is the median of the mass estimates derived from the fitting of the 300 Monte-Carlo perturbed SEDs of each class; and ${\mathit{M}}_{\mathrm{dust}}^{\mathrm{actual}}$ which is the dust mass corresponding to the unperturbed SED. In other words, ${\mathit{M}}_{\mathrm{dust}}^{\mathrm{measured}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}^{\mathrm{actual}}$ quantifies the deviation from the dust mass derived from an SED fit, to its actual value, as a function of the starlight intensity.

Figure B.2 shows that up to  ⟨ U ⟩ ^Std ≲ 6 and  ⟨ U ⟩ ^AC ≲ 9, the dust masses derived from SED fits are not significantly biased (i.e. the deviation is smaller than the error bars). However, for starlight intensities higher than these values, the SED fit tends to systematically underestimate the dust masses (i.e. overestimate the gas-to-dust mass ratio). The amplitude of this effect can go up to a factor of  ≃2, for high  ⟨ U ⟩ .

This is a demonstration of the effect invoked to explain the low gas-to-dust mass ratio, at high starlight intensities, in Fig. 14. Fortunately, this effect concerns only a regime containing a small fraction of the pixels. Moreover, these pixels are the less massive, since they correspond to hot and diffuse regions. Therefore, this bias does not have a significant impact on our global dust mass estimate.

Appendix C: Relevance of our starlight intensity distribution

C.1. The unnecessariness of adding a diffuse field component

We note that our approach of modelling IR SEDs with an empirical combination of starlight intensities, is common in the literature. In particular, Draine et al. (2007) modelled the SEDs of the SINGS galaxies with an extra, uniformly illuminated, component: $\begin{array}{ccc}\frac{\mathrm{d}{\mathit{M}}_{\mathrm{dust}}^{\mathrm{extra}}}{\mathrm{d}\mathit{U}}& \mathrm{=}& \mathit{\gamma }{\mathit{M}}_{\mathrm{dust}}\mathrm{\times }\frac{\mathrm{(}\mathit{\alpha }\mathrm{-}\mathrm{1}\mathrm{)}}{{\mathit{U}}_{\min }^{\mathrm{1}\mathrm{-}\mathit{\alpha }}\mathrm{-}\mathrm{(}{\mathit{U}}_{\min }\mathrm{+}\mathrm{\Delta }\mathit{U}{\mathrm{)}}^{\mathrm{1}\mathrm{-}\mathit{\alpha }}}{\mathit{U}}^{\mathrm{-}\mathit{\alpha }}\\ & \mathrm{+}& \mathrm{(}\mathrm{1}\mathrm{-}\mathit{\gamma }\mathrm{)}{\mathit{M}}_{\mathrm{dust}}\mathrm{\times }\mathit{\delta }\mathrm{(}\mathit{U}\mathrm{-}{\mathit{U}}_{\min }\mathrm{)}\mathit{,}\end{array}$(C.1)where γ is an extra parameter controlling the mass fraction of the U^−α component. This formulation implicitly assumes that there are no massive quantities of dust colder than the diffuse ISM. One of the advantages of this extra component is to avoid the dust mass to diverge, in the absence of submm constraints. Indeed, as demonstrated by Galametz et al. (2011), modelling the SEDs of galaxies without submm data, using Eq. (4), leads to gross errors on the dust mass. However, in our case, we have the valuable submm constraints provided by Herschel. The relevance of this extra component can be tested using the Monte-Carlo method, described in Sect. 3.4.2.

First, we tested the stability of this extra component. Figure C.1 shows the distribution of the Monte-Carlo values of γ, when fitting the integrated SED strip (R10), with Eq. (C.1) and the “standard” grain composition. This distribution is drastically bimodal. Most of the perturbed SEDs are fit either without the extra δ(U − U_min) components (γ ≃ 1) –  which is equivalent to our model – or with a sole uniformly illuminated SED (γ ≃ 0). In other words, the value of the parameter γ is totally uncertain.

Second, we tested the statistical relevance of adding this extra component. To quantify this aspect, we performed a F-test. We computed the statistic F_χ (Bevington & Robinson 2003): ${\mathit{F}}_{\mathit{\chi }}\mathrm{=}\frac{{\mathit{\chi }}_{\mathrm{Std}}^{\mathrm{2}}\mathrm{-}{\mathit{\chi }}_{\mathrm{extra}}^{\mathrm{2}}}{{\mathit{\chi }}_{\mathrm{extra}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathit{n}\mathrm{-}{\mathit{m}}_{\mathrm{extra}}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(C.2)where ${\mathit{\chi }}_{\mathrm{Std}}^{\mathrm{2}}$ is the χ² (Eq. (8)) of the fit of a given perturbed SED of the integrated strip with the “standard model”, and ${\mathit{\chi }}_{\mathrm{extra}}^{\mathrm{2}}$ is the corresponding value for Eq. (C.1); n = 10 is the number of wavebands, and m_extra = 7 is the number of free parameters for the model of Eq. (C.1) plus the stellar component of Eq. (7). F_χ can be seen as a measure of how much the additional term has improved the value of the reduced χ². Statistically, ${\mathit{F}}_{\mathit{\chi }}\mathrm{\simeq }\mathrm{0.05}{\mathrm{2}}_{-0.027}^{\mathrm{+}\mathrm{0.023}}\mathrm{\ll }\mathrm{0.67}$ is lower than the value of the F-distribution, with n − m_extra − 1 = 2 degrees of freedom, with a probability of exceeding F_χ of 50% (or P_F(0.67,1,n − m_extra − 1) = 50%; Bevington & Robinson 2003). In other words, adding the extra component of Eq. (C.1) does not significantly improve the χ². Our starlight intensity distribution (Eq. (4)) is therefore statistically more significant.

Finally, we tested the conservativeness of the dust mass estimate with Eq. (C.1). We have performed the fits of exactly the same perturbed SEDs that were used in Fig. C.1, with our “standard model” (Eq. (4)). The ratio of the mass of our model to the mass obtained with the extra component is ${\mathit{M}}_{\mathrm{dust}}^{\mathrm{Std}}\mathit{/}{\mathit{M}}_{\mathrm{dust}}^{\mathrm{extra}}\mathrm{\simeq }\mathrm{1.00}{\mathrm{0}}_{-0.045}^{\mathrm{+}\mathrm{0.010}}$. In other words, the two starlight intensity distributions of Eqs. (4) and (C.1) give statistically identical dust masses, within a few percents.

In summary, for the analysis performed in this paper, the starlight intensity distribution of Eq. (C.1) is not physically motivated, neither statistically stable nor relevant, and gives identical results to Eq. (4). Our model has an appropriate balance between free parameters and observational constraints. It is flexible enough, but it does not lead to overinterpreting the data.

C.2. Comparison of our model with the isothermal approximation

As justified in Sect. 3, our model accounts for a distribution of equilibrium dust temperatures within each pixel (Eq. (4)). This formalism allows us to fit the submm slope more accurately than with a single temperature. As a consequence, our model predicts more mass than a single black body fit having the same dust opacity. Since the single black body approach is still widely used, even in the Herschel era, we have performed a systematic comparison in order to quantify the biases of such an approach.

To perform our comparison, we have fitted the R4 (54 pc) map with a single modified black body having the same grain opacity as our “standard model”: κ_abs(160 μm) = 1.4m²   kg^-1, and β = 2 (Fig. A.1). For each pixel, the mass and the temperature are free to vary, but β is kept fixed. We constrain this model with the MIPS_70 μm, MIPS_160 μm, SPIRE_250 μm and SPIRE_500 μm fluxes, weighted the same way as our complete model (Eq. (8)), except that the weight of the MIPS_70 μm flux is divided by 100. The purpose of this trick is the following: (i) when the peak of the SED is well constrained by MIPS_160 μm, SPIRE_250 μm and SPIRE_350 μm (low temperatures), the MIPS_70 μm flux has almost no impact on the χ²; (ii) when the peak is not well constrained (high temperatures), the MIPS_70 μm becomes important and temperature divergence is avoided. The mass for each pixel is noted ${\mathit{M}}_{\mathrm{dust}}^{\mathrm{1}\mathrm{BB}\mathrm{(}\mathit{\beta }\mathrm{=}\mathrm{2}\mathrm{)}}$.

Figure C.2 shows the ratio between the dust masses obtained with this isothermal fit and with our “standard model”, as a function of the mass averaged starlight intensity  ⟨ U ⟩ . It shows that for  ⟨ U ⟩  ≳ 1, the mass ratio is roughly constant between 60 and 80%, but for  ⟨ U ⟩  ≲ 1, the ratio drops down to less than 30%. This simulation demonstrates that an isothermal fit will induce a bias in dust masses with average temperature. The regions with high starlight intensities (⟨ U ⟩  ≳ 1) are generally diffuse regions and their SEDs are almost isothermal. This explains why the mass ratio does not vary much with starlight intensity in this range. The average value of the ratio is around 70% in these regions, which is reasonable considering the crudeness of the approximation. On the contrary, regions with low starlight intensities (⟨ U ⟩  ≲ 1) are generally dense regions. Since the mean free path of photons is much shorter in these regions, there is a significant mix of cold and hot components within each pixel. The isothermal approximation is not valid anymore in this range of  ⟨ U ⟩ , and the mass obtained with a single black body fit is biased. Namely, it drops by a factor of  ≃2 compared to its diffuse ISM value, down to  ≃30% of the value of the “standard model”. We obtain a similar trend comparing the “AC model” with a single black body fit having κ_abs(160 μm) = 1.6m²   kg^-1, and β = 1.7, except that the values of  ⟨ U ⟩  are systematically shifted by a factor of  ≃2 (as demonstrated in Fig. 5).

Finally, we note that the increase of dust mass at low starlight intensities compared to the isothermal approximation is unlikely a bias induced by the submm excess extending down to the SPIRE_350 μm band. Indeed, this excess is negligible in regions with low  ⟨ U ⟩  (Sect. 4.4). It is prominent only in regions with high  ⟨ U ⟩ , where the trend of Fig. C.2 is not significant.

Appendix D: Inconsistency of the very cold dust hypothesis

As discussed by Galliano et al. (2003, 2005) and Galametz et al. (2009, 2010), the SPIRE_500 μm excess could be attributed to very cold dust (VCD; ${\mathit{T}}_{\mathrm{eq}}^{\mathrm{VCD}}\lesssim \mathrm{10}\mathrm{K}$). Simply considering the shape of the SED, this explanation is not unlikely. Such a component would produce a change in the submm slope of the SED. Indeed, warm and cold dust are heated by the transmitted stellar light, in different environments, at different optical depths. Their temperature distribution is therefore continuous. On the other hand, in order to reach very cold temperatures (T_dust ≲ 10K), the dust has to be fully shielded from stellar radiation. In these conditions, the dominant heating sources are the collisions and the IR radiation (galaxy and CMB) which do not depend on the optical depth. It could reach temperatures lower than T_dust ≃ 5K (Galliano et al. 2003) only with difficulty. Very cold dust therefore corresponds to a change of regime in the dust heating. If it exists, its temperature should not significantly vary with the optical depth in the cloud. Consequently its emission would be seen as a roughly isothermal component at 5K ≲ T_dust ≲ 10K. With sufficient mass, it would produce an excess emission, with a change of slope, at wavelengths 300 μm ≲ λ ≲ 1mm.

Until now, we did not have the spatial resolution to explore this hypothesis. In this section, we do an order of magnitude estimate to test the likeliness of the VCD model. What follows is an update of the discussion at the end of Galliano et al. (2003). This is a complement to Sect. 4.4.2.

Let’s first assume that the excess originates in equilibrium grains having the UV-to-IR cross-section of silicates: $\begin{array}{ccc}& & \mathit{V}\hspace{0.17em}\mathrm{bandopacity:}{\mathit{\kappa }}_{\mathrm{abs}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{V}}\mathrm{\right)}\mathrm{\simeq }\mathrm{130}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\mathrm{k}{\mathrm{g}}^{-1}\\ & & {\mathrm{SPIRE}}_{\mathrm{500} \mathit{\mu }\mathrm{m}\mathrm{)}}\mathrm{opacity:}{\mathit{\kappa }}_{\mathrm{abs}}\mathrm{\left(}{\mathrm{SPIRE}}_{\mathrm{500} \mathit{\mu }\mathrm{m}}\mathrm{\right)}\mathrm{\simeq }\mathrm{0.23}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}\mathrm{k}{\mathrm{g}}^{-1}\mathit{,}\end{array}$(D.1)and a submillimeter opacity index β^VCD = 1, and an equilibrium temperature of ${\mathit{T}}_{\mathrm{eq}}^{\mathrm{VCD}}\mathrm{\simeq }\mathrm{10}$ K. Those are the most optimistic values. If we can invalidate them, then the VCD hypothesis will be unrealistic for lower temperatures, steeper submm opacities, and larger UV cross-sections.

Although the relative excess is stronger in diffuse regions, we find an excess in almost every pixel of the strip (Fig. 16). Let’s assume that VCD lies in the core of very dense spherical clumps in the diffuse ISM. Then, the ISRF they have to be shielded from is the typical ISRF of the diffuse ISM. For the more realistic model (“AC”), this is basically  ⟨ U ⟩  ≳ 3. The optical depth, assuming a slab extinction, to shield this dust, and to allow it to reach ${\mathit{T}}_{\mathrm{eq}}^{\mathrm{VCD}}$ is ${\mathit{A}}_{\mathrm{V}}^{\mathrm{VCD}}\mathrm{\simeq }\mathrm{2.9}$. Since these clumps have to be small, and that we see the excess everywhere, it means that each clump is unresolved, even at R1. Therefore, the diameter of these clumps has to be lower than the pixel size of R1, which is l_pix(R1) ≃ 10 pc.

This maximum size translates into a minimum density of ${\mathit{n}}_{\mathrm{H}}^{\min }$, in order to reach the required optical depth of ${\mathit{A}}_{\mathrm{V}}^{\mathrm{VCD}}\mathrm{\simeq }\mathrm{2.9}$: ${\mathit{n}}_{\mathrm{H}}^{\min }\gtrsim \frac{\frac{{\mathit{A}}_{\mathrm{V}}^{\mathrm{VCD}}}{\mathrm{1.086}}\mathrm{\times }{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}}{{\mathit{m}}_{\mathrm{H}}\mathrm{\times }{\mathit{\kappa }}_{\mathrm{abs}}\mathrm{\left(}{\mathit{\lambda }}_{\mathit{V}}\mathrm{\right)}\mathrm{\times }\frac{{\mathit{l}}_{\mathrm{pix}}\mathrm{\left(}\mathit{R}\mathrm{1}\mathrm{\right)}}{\mathrm{2}}}\mathrm{\simeq }\mathrm{2.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{H}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}\mathit{,}$(D.2)m_Hbeing the mass of an H atom. Notice that, since we are performing only an order of magnitude estimate, we use the slab extinction for a sphere.

With  ⟨ U ⟩  ≃ 3, the temperature of the shielding dust, around the core, is ${\mathit{T}}_{\mathrm{dust}}^{\mathrm{shield}}\mathrm{\simeq }\mathrm{21}\mathrm{K}$. Thus, a typical excess of r₅₀₀ ≃ 15% corresponds to a mass ratio, between the VCD core and the shielding ISM dust around, of: $\frac{{\mathit{M}}_{\mathrm{dust}}^{\mathrm{VCD}}}{{\mathit{M}}_{\mathrm{dust}}^{\mathrm{shield}}}\mathrm{\simeq }\frac{{\mathit{r}}_{\mathrm{500}}}{\mathrm{1}\mathrm{-}{\mathit{r}}_{\mathrm{500}}}\frac{{\mathit{l}}_{\mathit{\nu }}^{\mathrm{VCD}}\mathrm{(}{\mathrm{SPIRE}}_{\mathrm{500} \mathit{\mu }\mathrm{m}\mathrm{)}}}{{\mathit{l}}_{\mathit{\nu }}^{\mathrm{shield}}\mathrm{\left(}{\mathrm{SPIRE}}_{\mathrm{500} \mathit{\mu }\mathrm{m}\mathrm{)}}\mathrm{\right)}}\mathrm{\simeq }\mathrm{0.5}\mathit{,}$(D.3)l_νbeing the specific monochromatic power of the components. Therefore, for this picture to be correct, there should be at least ${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{shield}}$ of shielding ISM dust surface density, such that: ${\mathrm{\Sigma }}_{\mathrm{dust}}^{\mathrm{shield}}\gtrsim \frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\left(\frac{{\mathit{l}}_{\mathrm{pix}}\mathrm{\left(}\mathit{R}\mathrm{1}\mathrm{\right)}}{\mathrm{2}}\right)}^{\mathrm{3}}\frac{{\mathit{n}}_{\mathrm{H}}^{\min }{\mathit{m}}_{\mathrm{H}}}{{\mathit{G}}_{\mathrm{dust}}^{\exp \mathit{.}}}\frac{\mathrm{1}}{{\mathit{l}}_{\mathrm{pix}}\mathrm{(}\mathit{R}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}}\mathrm{\simeq }\mathrm{10}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}\mathrm{p}{\mathrm{c}}^{-2}\mathit{.}$(D.4)This value is 2 orders of magnitude higher than the typical dust mass surface density, and one order of magnitude higher than the highest values (Fig. 13). Moreover, this argument is conservative, since the highest excesses are found in lowest surface density regions (Sect. 4.4).

In summary, we have shown that the SPIRE_500 μm excess could not be accounted for by very cold dust, since it would require a minimum mass of shielding dust, within each pixel, too large compared to the observed surface density. It is equivalent to note that, the excess being present in most pixels of R1, we would need at least 10⁵ clumps of cold dust. It is more efficient to hide very cold dust in a small number of clumps. For example, Galliano et al. (2003) estimated that this number could be less than a few hundreds in the Magellanic dwarf galaxy NGC1569.

---

Acknowledgments

We would like to thank Anthony Jones, Vincent Guillet, Laurent Verstraete, Bruce Draine and Karl Gordon for stimulating discussions on the matter of this paper. We thank the anonymous referee, for his report, and the editor, Malcolm Walmsley, for his careful reading of the paper. We acknowledge extensive use of the Levenberg-Marquardt χ² minimization routine written by Craig Markwardt. SPIRE has been developed by a consortium of institutes led by Cardiff Univ. (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA).

References

All Tables

All Figures