Dust Abundance Variations in the Magellanic Clouds: Probing the Lifecycle of Metals with All-Sky Surveys

TL;DR

This study uses all-sky infrared surveys to measure how dust-to-gas ratios vary in the Magellanic Clouds, revealing significant decreases from diffuse to dense regions, which impacts galaxy evolution models.

Contribution

It provides the first all-sky survey-based measurements of dust abundance variations in the Magellanic Clouds, overcoming previous calibration limitations.

Findings

G/D decreases by factors of 3 in the LMC and 7 in the SMC from diffuse to dense ISM.

The G/D variation aligns with elemental depletion and metal accretion models.

Results inform galaxy evolution simulations and dust-based gas mass estimates.

Abstract

Observations and modeling suggest that the dust abundance (gas-to-dust ratio, G/D) depends on (surface) density. The variations of the G/D provide constraints on the timescales for the different processes involved in the lifecycle of metals in galaxies. Recent G/D measurements based on Herschel data suggest a factor 5---10 decrease in the dust abundance between the dense and diffuse interstellar medium (ISM) in the Magellanic Clouds. However, the relative nature of the Herschel measurements precludes definitive conclusions on the magnitude of those variations. We investigate the variations of the dust abundance in the LMC and SMC using all-sky far-infrared surveys, which do not suffer from the limitations of Herschel on their zero-point calibration. We stack the dust spectral energy distribution (SED) at 100, 350, 550, and 850 microns from IRAS and Planck in intervals of gas surface density, model the stacked SEDs to derive the dust surface density, and constrain the relation between G/D and gas surface density in the range 10---100 \Msu pc$^{-2}$ on $\sim$ 80 pc scales. We find that G/D decreases by factors of 3 (from 1500 to 500) in the LMC and 7 (from 1.5$\times 10^4$ to 2000) in the SMC between the diffuse and dense ISM. The surface density dependence of G/D is consistent with elemental depletions and with simple modeling of the accretion of gas-phase metals onto dust grains. This result has important implications for the sub-grid modeling of galaxy evolution, and for the calibration of dust-based gas mass estimates, both locally and at high-redshift.

Tables

Table 1. Table 1 Coefficients of the fits to the relation between MW H  i column density and cirrus dust emission

	LMC			SMC
	Linear	Quadratic		Linear	Quadratic
Band	$a_1$	$a_1$	$a_2$	$a_1$	$a_1$	$a_2$
100	0.681$\pm$8.9$\times {10}^{-4}$	0.642$\pm$2.8$\times {10}^{-3}$	0.00706$\pm$4.7$\times {10}^{-4}$	0.404$\pm$1.8$\times {10}^{-3}$	0.145$\pm$9.0$\times {10}^{-3}$	0.0802$\pm$2.7$\times {10}^{-3}$
350	0.436$\pm$5.5$\times {10}^{-4}$	0.362$\pm$1.6$\times {10}^{-3}$	0.0138$\pm$2.9$\times {10}^{-4}$	0.383$\pm$1.7$\times {10}^{-3}$	0.265$\pm$8.6$\times {10}^{-3}$	0.0357$\pm$2.6$\times {10}^{-3}$
550	0.149$\pm$1.9$\times {10}^{-4}$	0.129$\pm$5.6$\times {10}^{-4}$	0.00370$\pm$9.9$\times {10}^{-5}$	0.129$\pm$5.7$\times {10}^{-4}$	0.0897$\pm$2.9$\times {10}^{-3}$	0.0120$\pm$8.6$\times {10}^{-4}$
850	0.0693$\pm$8.8$\times {10}^{-5}$	0.0932$\pm$2.9$\times {10}^{-4}$	-0.00400$\pm$4.7$\times {10}^{-5}$	0.0769$\pm$3.4$\times {10}^{-4}$	0.102$\pm$1.7$\times {10}^{-3}$	-0.00726$\pm$4.8$\times {10}^{-4}$

Table 2. Table 2 SEDs, gas and depletion measurements used to calibrate κ 160 subscript 𝜅 160 \kappa_{160} in the LMC and SMC

	$S_{100}$	$S_{350}$	$S_{550}$	$S_{850}$	${\Sigma }_g$	$F_{\ast }$	${\left(G/D\right)}_{\mathrm{dep}}$	${\kappa }_{160}$
	MJy/sr	MJy/sr	MJy/sr	MJy/sr	$M_{\odot }$ pc-2			cm2 g-1
LMC	22.4$\pm$16.5	8.2$\pm$4.1	2.9$\pm$1.4	0.87$\pm$0.42	22	0.42	328	6.1
SMC	13.6$\pm$1.05	5.95$\pm$0.43	2.37$\pm$0.16	0.80$\pm$0.051	65	0.80	800	1.65

Table 3. Table 3 Total dust and gas masses in the LMC and SMC

	LMC			SMC
	Gas	Dust	G/D	Gas	Dust	G/D
	$M_{\odot }$	$M_{\odot }$		$M_{\odot }$	$M_{\odot }$
Quadratic fit to MW cirrus emission	(4.0$\pm$0.04)$\times {10}^8$	(5.1$\pm$0.47)$\times {10}^5$	785 $\pm$ 72	(4.7$\pm$0.04)$\times {10}^8$	(7.3$\pm$0.77)$\times {10}^4$	6380 $\pm$ 680
Linear fit to MW cirrus emission	(4.0$\pm$0.04)$\times {10}^8$	(5.3$\pm$0.5)$\times {10}^5$	758 $\pm$ 68	(4.8$\pm$0.04)$\times {10}^8$	(6.7$\pm$0.7)$\times {10}^4$	7070 $\pm$ 770

Equations

$$S_{\nu }(\lambda )=(2.0891\times 10^4){\kappa }_0{\left(\frac{{\lambda }_0}{\lambda }\right)}^{\beta }{\Sigma }_d{B}_{\nu }(\lambda ,T_d)$$

$$S_{\mathrm{band}}^{\mathrm{mod}}=\frac{\int R_E(\nu )S_{\nu }(\nu )d\nu }{\int R_E(\nu )\left(\frac{{\nu }_{\mathrm{eff}}}{\nu }\right)d\nu }$$

$$P({\log }_{10}{\Sigma }_d,T_d,\beta |S_{\mathrm{band}}^{\mathrm{mean}})=\frac{1}{Q}\exp \left(-\frac{1}{2}{\chi }^2\right)$$

$${\chi }^2=\left(S_{\mathrm{band}}^{\mathrm{mod}}-S_{\mathrm{band}}^{\mathrm{mean}}\right){\mathbb{C}}^{-1}\left(S_{\mathrm{band}}^{\mathrm{mod}}-S_{\mathrm{band}}^{\mathrm{mean}}\right)$$

$$Q=(2\pi {)}^n\det |\mathbb{C}|$$

$${\mathbb{C}}_{ij}^{\mathrm{bkg}}=\frac{{\sum }_{k=0}^{N_r-1}\left(S_i^k-<S_i>\right)\left(S_j^k-<S_j>\right)}{N_{\mathrm{stack}}{N}_r}$$

$${\footnotesize\bf{\mathbb{C}^{\mathrm{bkg}}}=\frac{1}{N_{\mathrm{stack}}}\left[\begin{array}[]{cccc}2.51&0.38&0.14&0.040\\ 0.38&0.12&0.046&0.013\\ 0.14&0.046&0.018&0.0053\\ 0.040&0.013&0.0076&0.0018\\ \end{array}\right]}$$

$${\footnotesize\bf{\mathbb{C}^{\mathrm{bkg}}}=\frac{1}{N_{\mathrm{stack}}}\left[\begin{array}[]{cccc}0.088&0.035&0.013&0.0035\\ 0.035&0.020&0.0080&0.0025\\ 0.013&0.0080&0.0035&0.0011\\ 0.0035&0.0025&0.0011&0.00042\\ \end{array}\right]}$$

$${\footnotesize\bf{A}=\left[\begin{array}[]{cccc}1.036\times 10^{-2}&0&0&0\\ 0&4.096\times 10^{-3}&4.000\times 10^{-4}&0\\ 0&4.000\times 10^{-4}&3.721\times 10^{-3}&0\\ 0&0&0&6.0834\times 10^{-5}\\ \end{array}\right]}$$

$$P({\log }_{10}{\Sigma }_d|S_{\mathrm{band}}^{\mathrm{mean}})={\int }_{T_d,\beta }P({\log }_{10}{\Sigma }_d,T_d,\beta |S_{\mathrm{band}}^{\mathrm{mean}})$$

$${\left({\log }_{10}{\Sigma }_d\right)}^{\mathrm{best}}={\int }_{{\Sigma }_d}P({\log }_{10}{\Sigma }_d|S_{\mathrm{band}}^{\mathrm{mean}}){\Sigma }_{d}d{\log }_{10}{\Sigma }_d$$

$$\begin{gathered} \sigma {\left({\log }_{10}{\Sigma }_d\right)}^2={\int }_{{\Sigma }_d}P({\log }_{10}{\Sigma }_d|S_{\mathrm{band}}^{\mathrm{mean}})\times \\ {\left({\log }_{10}{\Sigma }_d-{\left({\log }_{10}{\Sigma }_d\right)}^{\mathrm{best}}\right)}^{2}d{\log }_{10}{\Sigma }_d \end{gathered}$$

$${\Sigma }_d^{\mathrm{fit}}=2.5\times 10^{-4}{\Sigma }_g^{1.5}$$

$$\begin{gathered} {\Sigma }_d^{\mathrm{fit}}=4.87\times 10^{-4}+4.19\times 10^{-5}{\Sigma }_g \\ +8.18\times 10^{-7}{\Sigma }_g^2+2.06\times 10^{-8}{\Sigma }_g^3 \end{gathered}$$

$$\begin{gathered} \left(\frac{{\tau }_{acc}}{yr}\right)=2\times 10^7\times \\ \left(\frac{<a>}{0.1\mu m}\right){\left(\frac{n_H}{100{\mathrm{cm}}^{-3}}\right)}^{-1}{\left(\frac{T_d}{50K}\right)}^{-\frac{1}{2}}{\left(\frac{Z}{0.02}\right)}^{-1} \end{gathered}$$

$$f_d=f_0\frac{\exp \left(\frac{t}{{\tau }_{acc}}\right)}{1-f_0+f_0\exp \left(\frac{t}{{\tau }_{acc}}\right)}$$

$$n_H=n_0\frac{{\Sigma }_g}{{\Sigma }_0}$$

$$\left(\frac{G}{D}\right)={\left(\frac{G}{D}\right)}_0\frac{f_0}{f_d}$$

Full text

1. Introduction

Dust grains absorb stellar light in the UV-optical and re-emit it in the FIR, which represents 30—50% of the output of a galaxy. Therefore, our ability to interpret observations of galaxies and trace their stellar, dust, and gas content over cosmic time over the entire spectral range critically relies on our understanding of how the dust abundance and properties vary with environment; this in turn requires us to understand the processes responsible for dust formation, destruction, and transport and their associated timescales. Addressing the question of dust growth and destruction requires large-scale maps of the dust (in the far-infrared, FIR) and gas (HI and CO) content of galaxies resolved on kpc-scales and over a range of metallicities. Subsequent gas-to-dust ratio and dust property measurements can then be compared to phenomenological, resolved models of the chemical and dust evolution of entire galaxies. These models have advanced to the point of needing more detailed, environment-dependent observations of dust abundances and properties. Since variations in the dust abundance and properties cannot be observationally characterized in unresolved galaxies, we must therefore characterize and understand the micro-astrophysical processes responsible for dust evolution in the ISM in nearby galaxies, and apply our knowledge to the more distant universe.

We have been investigating dust evolution in the Magellanic Clouds (LMC and SMC) by measuring their dust and gas contents using over 1400 hours of space observatory time with Spitzer (Surveying the Agents of Galaxy Evolution, SAGE survey Meixner et al., 2006; Gordon et al., 2011) and Herschel (HERschel Inventory of The Agents of Galaxy Evolution, HERITAGE survey Meixner et al., 2013a), in addition to ground-based measurements of CO (Mizuno et al., 2001b, a; Wong et al., 2011), HI (Stanimirovic et al., 1999; Kim et al., 2003a), and H$\alpha$ (Gaustad et al., 2001). Thanks to these observations, we now have dust surface density and temperature maps in the LMC and SMC (Gordon et al., 2014). We have developed sophisticated tools to extract accurate dust parameters (surface density, temperature etc.) from the images using probabilistic techniques (see Gordon et al., 2014, for details). In particular, systematics that can potentially affect dust and gas measurements in the low signal-to-noise, low surface density regime were carefully addressed and removed. For instance, covariance between different bands was accounted for; the dust emissivity was self-calibrated on Milky Way observations; and a foreground subtraction was done to remove the structured emission due to Milky Way (MW) dust (cirrus) emission. This step was particularly important for the SMC, where foreground Milky Way structures with similar surface brightness to those in the SMC were removed.

In Roman-Duval et al. (2014), we examined the resolved relation between dust and gas surface densities in the LMC and SMC. We discovered that 1) the hydrogen-to-dust ratio (H/D) in the diffuse atomic ISM is 380${}^{+250}_{-130}$ in the LMC and 1200${}^{+1600}_{-420}$ in the SMC; 2) the dust abundance (dust-to-hydrogen ratio, D/H) and FIR emissivity increase by a factor 2—3 between the diffuse atomic and the dense molecular ISM; 3) there is a gas pedestal of 10 $M_{\odot}$ pc*-2* (LMC) and 40 $M_{\odot}$ pc*-2* (SMC) in the relation between dust and hydrogen gas, implying the presence of a dust-poor ISM component at surface densities barely or not detected by Herschel (below the gas pedestal) located in the outskirts of the galaxies and holes in the ISM. This dust-poor ISM component would have a D/H 5—10 times lower than in the FIR detected regions and correspond to the junction between the ISM and the circum-galactic medium (CGM). This hypothesis explaining the non-zero intercept of the dust-gas relation is supported by a stacking analysis of the pixels of the FIR maps where the brightness is below the sensitivity cut and a dust surface density cannot be derived in Gordon et al. (2014). The dust mass corresponding to this FIR-faint component is (5.9$\pm$3.6) $\times$ 104 $M_{\odot}$ in the LMC and (1.6$\pm$1.3) $\times$ 104 $M_{\odot}$ in the SMC. The corresponding H i masses are 1.6 $\times$ 108 $M_{\odot}$ and 2.0 $\times$ 108 $M_{\odot}$ in the LMC and SMC, leading to global H/D of $\sim$2800 in the LMC and $\sim$12,000 in the SMC, or 5—10 times higher than in the diffuse atomic ISM detected in the FIR.

While this tantalizing result could have important implications for the chemical evolution of the LMC and SMC, and galaxies in general, its robustness is limited by the inherently relative nature of the Herschel (and Spitzer) observations. Indeed, in order to remove 1/f noise (striping) and instrumental drifts from individual scans, one has to make assumptions about the level and structure of the emission in the periphery of the observed galaxies. For observations of galaxies relatively small on the sky, observations usually include a suitable background area far from the center of the galaxy where zero-emission can reasonably be assumed. For close-by galaxies however, this turns out to be problematic since the observations usually do not extend far enough. As a result, although there is still a significant level of astrophysical signal at the periphery of the images, this low-level emission is zeroed out in the maps. This is illustrated in Figure 1, which shows a difference image between the 350 $\mu$m bands of Planck and Herschel/SPIRE at 5’ resolution. To circumvent this issue in Roman-Duval et al. (2014), we performed a background subtraction on the H i 21 cm images to put them on the same pedestal as the dust maps. While this provides a reasonable data set to evaluate variations in the dust abundance, an absolute measurement of the gas-to-dust ratio is not possible. Additionally, issues remain due to the different background assumptions between the scans and cross-scans in the Herschel images.

In this paper, we use Planck and IRAS all-sky observations, which are not affected by this limitation on the zero-point calibration and cover the peak of the dust SED at 100, 350, 550, and 850 $\mu$m, to quantify the spatial variations of the dust abundance, from the very diffuse ISM in the outskirts of the LMC and SMC, to the dense ISM where CO emission is detected. We separate the ISM in intervals of HI column density, and stack the dust emission in each gas surface density bin. We then fit a dust model to the resulting stacked SED to determine the mean dust surface density, temperature, and spectral emissivity index in each bin. We then compare the resulting trend in D/G vs $\Sigma_{g}$ to resolved chemical evolution models. In addition to improvements on the zero-points, this stacking technique allows us to probe the dust content of lower surface density regions compared to Roman-Duval et al. (2014), and to measure the variations of G/D as a function of surface density.

The paper is organized as follows. In Section 2, we describe the data sets used in this analysis. In Section 3, we describe the methodology used to subtract the foreground Milky Way cirrus emission, the stacking of the dust emission in intervals of gas surface density, and the dust modeling. In Section 4, we present the main results of this paper, specifically the variations of the dust surface density, gas-to-dust ratio, and other dust parameters with gas surface density. Finally, we discuss the possible origins of the observed trends and their implications in Section 5.

2. Observations

The observations used in this analysis are described in this section, and include FIR dust emission observed in IRAS and PLANCK, 21 cm emission from H i gas observed with both single dish and interferometry, and CO emission from molecular gas observed with various ground-based millimeter facilities. Figure 2 shows the IRAS 100 $\mu$m (blue), Planck 350 $\mu$m (green), and H i surface density (red) maps as three-color images.

2.1. Far-infrared dust observations

This work utilizes all-sky observations with IRAS 100 $\mu$m and PLANCK 350, 550, and 850 $\mu$m. We used the ISSA (IRAS Sky Survey Atlas) version of the IRAS observations, which were downloaded from the Infra-Red Science Archive (IRSA). The Improved reprocessing of the IRAS Survey (IRIS) is affected by a $\sim$ 30% gain offset when compared to COBE/DIRBE and Herschel in the area of the Magellanic Clouds (Meixner et al., 2013b), and we hence chose to use the older ISSA version.

The Planck maps at 350, 550, and 850 $\mu$m (857 GHz, 545 GHz, 353 Ghz) from the 2nd release (Planck Collaboration et al., 2016a) were downloaded in HEALPIX format from the Planck archive. They were convolved to 5’ resolution, re-projected in WCS format, and converted to MJy/sr. Cosmic microwave background (CMB) emission and fluctuations were subtracted using the map from Bobin et al. (2016), which is the only CMB map not impainted in regions like the Magellanic Clouds and does not show significant residuals of foreground emission in bright regions.

Cosmic infrared background (CIB) emission was subtracted using mean intensity levels of 0.11 MJy/sr, 0.64 MJy/sr, 0.35 MJy/sr, and 0.13MJy/sr at 100 $\mu$m, 350 $\mu$m, 550 $\mu$m, and 850 $\mu$m, respectively (Planck Collaboration et al., 2016c, a).

2.2. Atomic gas

2.2.1 Observations

The atomic gas surface density is derived from 21 cm line emission maps with 60” and 98” resolution in the LMC and SMC. The HI maps combine observations carried out at the Australian Telescope Compact Array (ATCA) and at the Parkes 64 m radio Telescope. The combination of the single dish and interferometric observations are described in Kim et al. (2003b, LMC) and Stanimirovic et al. (1999, SMC). The H i column density maps are converted from the 21 cm emission under the assumption of an optically thin line (e.g. Bernard & et al., 2008). This assumption may not be correct in regions of high H i column density, where we may systematically underestimate the H i content (Dickey et al., 2000a; Fukui et al., 2014).

2.2.2 Conversion between 21 cm emission and H i surface density

The H i column density in the ATCA+Parkes survey, $N($H i$)$ in units of cm*-2*, is converted to a surface density, $\Sigma$(H i), via $\Sigma$(H i) $=$ $0.8\times 10^{-20}$ $N($H i$)$, where $\Sigma$(H i) is in $M_{\odot}$ pc*-2*. This conversion does not include the mean molecular weight of Helium. The atomic gas surface density, including Helium, is given by $\Sigma_{g,\mathrm{atomic}}$ $=$ 1.36 $\Sigma$(H i).

2.2.3 H i maps of Milky Way Foreground

To estimate MW cirrus FIR emission, we use the atomic gas surface density images in the MW velocity range. For the SMC, these are obtained from the ATCA+Parkes cubes (Stanimirovic et al., 1999), integrated over the MW velocity range. For the LMC, the MW velocities in the ATCA+Parkes data (Kim et al., 2003b) were never reduced, and we make use of the Parkes only cube instead (Staveley-Smith et al., 2003), integrated over MW velocities.

2.2.4 Stray-light correction

The ATCA+Parkes 21 cm data are not corrected for stray-light, which may contribute a low-level of extended emission. To estimate the contribution of stray-light to the H i surface density measured in the ATCA+Parkes data, we use the corrected GASS survey (16’ resolution McClure-Griffiths et al., 2009; Kalberla & Haud, 2015), which was downloaded from https://www.astro.uni-bonn.de/hisurvey/gass/. While the GASS survey has too coarse a resolution (16’ corresponding to 230 pc in the LMC and 290 pc in the SMC) to be used to derive the G/D trends, it can be used to estimate the contribution of stray-light in the ATCA+Parkes maps. We convolve the ATCA+Parkes maps to GASS resolution using Gaussian kernels and resample the convolved ATCA+Parkes and GASS maps on the same astrometric grid, as described in Section 2.4. We examine the collapsed HI 21 cm spectra in Y-log scale and determined the atomic gas velocity range of the LMC and SMC to be 150—600 km s*-1* and 55-270 km s*-1*, respectively. We subsequently extract Milky Way (MW), LMC, and SMC HI 21 cm spectral cubes and integrated intensity images for each galaxy. The cubes were generated in units of K, and so the integrated intensity images are in K km s*-1*. We convert the intensity to a column density as in Equation 2 of Bernard & et al. (2008) via $N($H i$)$ $=$ 1.82$\times 10^{18}$ $W($H i$)$, where $N($H i$)$ is the H i column density in cm*-2* and $W($H i$)$ is the 21 cm integrated intensity in K km s*-1*. We compute the difference between the H i column densities derived from the 16’ ATCA+Parkes and GASS maps, which corresponds to the contribution of stray-light on scales larger than 16’ in the ATCA+Parkes data. The correction, shown in Figure 3, amounts to an equivalent H i surface density $<$ 10 $M_{\odot}$ pc*-2* in the LMC and $<$ 4 $M_{\odot}$ pc*-2* in the SMC. Finally, we resample the stray-light correction to the astrometric grid of the ATCA+Parkes, which matches that of the other maps (Planck, IRAS, CO), and subtract it from the ATCA+Parkes observations.

2.3. Molecular gas traced by CO

The molecular gas component is partially traced by its 12COJ = 1-0 emission in the NANTEN survey of the LMC and SMC (Mizuno et al., 2001a) at 2.6*′* resolution. The NANTEN surveys of the Magellanic Clouds have complete coverage and should capture most of the molecular gas content traceable via CO emission. At this resolution, we assume a CO-to-molecular gas conversion factor $\alpha_{\mathrm{CO}}$ $=$ 6.4 $M_{\odot}$ pc*-2* (K km s${}^{-1})^{-1}$ ($\sim$1.5$\times$ Galactic) in the LMC, and $\alpha_{\mathrm{CO}}$ $=$ 21 $M_{\odot}$ pc*-2* (K km s${}^{-1})^{-1}$ ($\sim$5$\times$ Galactic) in the SMC (Bolatto et al., 2013). This conversion factor includes the mean molecular weight of Helium. The molecular gas surface density is then $\Sigma_{g,\mathrm{mol}}$ $=$ $\alpha_{\mathrm{CO}}$ $I_{\mathrm{CO}}$, where $I_{\mathrm{CO}}$ is the 12CO integrated intensity (in K km s*-1*).

2.4. Convolution to common resolution and resampling

All maps were convolved to the 5’ limiting resolution of the Planck observations, corresponding to 70 pc in the LMC and 90 pc in the SMC . We used a gaussian kernel, of FWHM equal to the quadratic difference between the final and original resolutions, to perform the convolution. The maps were then resampled onto a common astrometric grid of pixel size 5’ to yield independent pixels using the IDL routine HASTROM.

3. Method

The methodology used to derive the dust abundance (G/D or D/G) as a function of gas surface density is as follows. First, MW cirrus emission is estimated and removed from the FIR maps. The level of FIR emission by foreground MW cirrus can be comparable to faint emission in the periphery of the LMC or SMC, and so subtracting the foreground cirrus is particularly important for our purpose of measuring the dust abundance in very diffuse regions of the LMC and SMC. Second, the FIR dust emission is summed (or ”stacked”) over pixels in a given interval of gas surface density. A modified black body model is fit to the stacked SED to obtain the dust surface density, dust temperature, and spectral emissivity index for that gas surface density. The dust abundance (G/D) is then computed as the ratio of the gas to the dust surface density in each interval.

3.1. Benefits of the stacking analysis on all-sky surveys

In Roman-Duval et al. (2014), we used the dust surface density maps derived from Herschel data in Gordon et al. (2014) convolved to 1’ (15 pc) in the LMC and 98” (45 pc) in the SMC to examine the gas-to-dust ratio in the atomic and molecular ISM. Roman-Duval et al. (2014) performed a linear fit to the relation between dust and gas surface densities in the regime where 1) gas is completely atomic and 2) molecular gas is traced by CO, and thus obtained a single G/D value for each regime. Figure 4 shows the distribution of H i surface densities in the analysis by Roman-Duval et al. (2014). The linear fits to the dust-gas relation are dominated by H i surface densities $\Sigma$(H i) $=$ 10—20 $M_{\odot}$ pc*-2* ($\Sigma_{g}$ $=$ 14—27 $M_{\odot}$ pc*-2*) in the LMC and $\Sigma$(H i) $=$ 40—60 $M_{\odot}$ pc*-2* ($\Sigma_{g}$ $=$ 55—80 $M_{\odot}$ pc*-2*) in the SMC. The distribution of H i surface densities used in this analysis based on stacking of the Planck and IRAS data is also shown in Figure 4, and extends down to much lower H i surface density in the SMC over a much larger area (as shown by the number of included pixels). Similarly in the LMC, a substantially larger area in diffuse regions is included in this new analysis. The additional coverage is illustrated in Figure 2, which shows a 3-color image of the LMC and SMC with the area included in each analysis.

In addition to improvements on the zero-point calibration (see Section 1), the stacking analysis performed in this paper is thus focused on the variations of G/D as a function of surface density, down to surface densities substantially lower than Roman-Duval et al. (2014), and on scales larger than in Roman-Duval et al. (2014) (75 pc/90 pc in the LMC/SMC here vs 15 pc/45 pc in Roman-Duval et al., 2014). Ultimately, we were able to measure the variations of G/D with surface density between $\Sigma_{g}$ $\sim$ 8 $M_{\odot}$ pc*-2* and 100 $M_{\odot}$ pc*-2* on 75 pc/90 pc scales, instead of a bimodal estimation of G/D in the atomic and CO-bright ISM in Roman-Duval et al. (2014).

Lastly, we have considerably improved the subtraction of foreground MW cirrus emission compared to Gordon et al. (2014). In Gordon et al. (2014), the MW cirrus subtraction was done by converting the MW H i column density toward the LMC and SMC to FIR emission assuming conversion factors modified from Compiègne et al. (2011), derived from the average diffuse MW SED. This resulted in significant over-subtraction in the outskirts of the LMC and SMC, likely due to the fact that the color of the MW cirrus varies with position. Here, we empirically derive the conversion factors specifically for each region (LMC and SMC), which allows us to reduce the systematic uncertainties in our FIR maps and accurately estimate the diffuse dust emission from the LMC and SMC.

3.2. Milky Way Cirrus Subtraction

MW dust cirrus emission occurs in the foreground of the LMC and SMC dust emission at a level of a few MJy/sr at the peak of the SED (in the range 100-350 $\mu$m). We estimate the detailed structure of the MW dust emission in the IRAS and Planck bands using the surface density map of H i in the MW, $\Sigma_{\mathrm{MW}}$(H i), in the direction of the LMC (Staveley-Smith et al., 2003) and SMC (Stanimirovic et al., 1999). The MW H i surface density was obtained by integrating the H i cubes over MW velocities. The cirrus dust emission, $S_{\nu}^{\mathrm{MW}}$ in MJy/sr, in each band was obtained by applying either a linear ($S_{\nu}^{\mathrm{MW}}=a_{1}\Sigma_{\mathrm{MW}}$(H i)) or quadratic ($S_{\nu}^{\mathrm{MW}}=a_{1}\Sigma$(H i) $+$ $a_{2}\Sigma_{\mathrm{MW}}$(H i)2) function to the MW H i surface density maps, $\Sigma_{\mathrm{MW}}$(H i). We constrained the linear and quadratic coefficients ($a_{1}$ and $a_{2}$) by examining the relation between $S_{\nu}^{\mathrm{MW}}$ and $\Sigma_{\mathrm{MW}}$(H i) in ”cirrus calibration regions”, where no H i emission at LMC or SMC velocities is detected to within the sensitivity limits of our maps, corresponding to pixels with H i surface densities below 4 $M_{\odot}$ pc*-2*. The contours corresponding to this threshold are shown in Figure 5. In those regions (outside the contours), the paucity of gas originating in the LMC/SMC guarantees the absence of dust emission from those galaxies, particularly given their low metallicty. Therefore, dust emission in those regions originates solely from Galactic foreground, and $S_{\nu}$ $=$ $S_{\nu}^{\mathrm{MW}}$, where $S_{\nu}$ is the emission in the maps. Figure 6 shows the relation between MW H i surface density and cirrus dust emission in each band and toward each galaxy, as well as the linear and quadratic fits to those relations. The coefficients of the fits are given in Table 1.

We note that the cirrus emission exhibits noticeable differences in color between the sight-lines toward the LMC and SMC, particularly in colors involving the 100 $\mu$m band. Hence, an attempt at fitting the relation between MW H i and dust emission simultaneously for both sight-lines failed, because it resulted in an over-subtraction in the IRAS 100 $\mu$m band in the SMC. Hence, we estimate the cirrus emission in the LMC and SMC sight-lines separately.

Figure 6 shows that a quadratic function provides a better description of the relation between MW H i surface density and cirrus dust emission. In the following, we therefore use the maps obtained from MW cirrus subtraction based on a quadratic fit as our fiducial case. This result is not unexpected. As the gas surface density increases, the temperature and spectral emissivity index of the dust vary, resulting in a non-linear relation between gas surface density and FIR emission.

3.3. Stacking of the dust emission

Once the foreground MW emission is removed, we proceed with the stacking analysis of the dust emission. The stacking analysis allows us to beat the noise in diffuse regions with faint FIR emission and measure dust emission at surface densities much lower than previously observed in Gordon et al. (2014) and Roman-Duval et al. (2014). Since dust and gas are mixed in the ISM, it is expected that the dust emission correlates best with gas surface density. We therefore stack the dust emission in bins of gas surface density.

The pixels of the dust maps are separated into bins of gas surface density. Each bin of gas surface density is 0.05 dex wide. The number of pixels in each bin, shown in Figure 7, varies between 5 and 2000. The decreasing number of pixels in bins of increasing surface density is due to the hierarchical (fractal) nature of the ISM. The diffuse, low surface density ISM occupies a much bigger volume that its dense counterpart. In order to sample the high range of the surface density, one has to use bins of equal surface density width, rather than the statistically advantageous equal number of pixels in each bin.

For each gas surface density bin, we compute the mean, standard deviation, and error on the mean of the dust emission in each band, $S_{\mathrm{band}}^{\mathrm{mean}}$, $S_{\mathrm{band}}^{\mathrm{stdv}}$, and $S_{\mathrm{band}}^{\mathrm{err}}$ $=$ $S_{\mathrm{band}}^{\mathrm{stdv}}$/$\sqrt{N_{\mathrm{stack}}}$, from all the $N_{\mathrm{stack}}$ pixels in each gas surface density interval. We have verified that the mean gas surface density in each bin equals the center of the surface density bin to within 5%. The stacked dust emission in each band is plotted as a function of gas surface density in Figure 8.

3.4. Modeling of the dust emission

We then fit a modified black body model (SMBB model in Gordon et al., 2014) to the stacked dust emission SED, $S_{\mathrm{band}}^{\mathrm{mean}}$, at 100, 350, 550, 850 $\mu$m in each bin of gas surface density. The modified black body function takes the form:

[TABLE]

where $S_{\nu}$ is the surface brightness, $B_{\nu}(\lambda,T_{d})$ is the Planck function (and $S_{\nu}$ and $B_{\nu}(\lambda,T_{d})$ have the same dimensions and units), $\kappa_{0}$ is the dust emissivity of the big grains at wavelength $\lambda_{0}$ in cm2 g*-1*, and $\Sigma_{d}$ is the dust surface density in $M_{\odot}$ pc*-2*. $\Sigma_{d}$, $T_{d}$, and $\beta$ are left as free parameters. We note that we also performed the dust SED modeling with the model of a black body modified by a broken-emissivity law (BEMBB model described by Equations 7, 9, and 10 in Gordon et al., 2014). While the dust parameters we derived (surface density, abundance, temperature, spectral emissivity index) are consistent between the two modeling approaches, uncertainties on the derived parameters were significantly larger for the BEMBB model due to the numbers of bands (4 in this study) being lower than the number of parameters (5 for the BEMBB model). Given the lack of scientific benefit in using the BEMBB model (same results, larger error bars) and, as we will see, the absence of significant residuals with the SMBB model, we only use the SMBB model described in Equation 1 in the rest of this paper.

We calibrate $\kappa_{0}$ at $\lambda_{0}$ $=$ 160 $\mu$m using the Milky Way diffuse dust SED from Compiègne et al. (2011) and the knowledge of the dust-to-gas ratio toward this sight-lines from depletions (Jenkins, 2009). The self-calibration of $\kappa_{0}$ is explained in details in Section 3.5, and leads to $\kappa_{0}$ $=$ 12.4 cm2 g*-1* at $\lambda_{0}$ $=$ 160 $\mu$m, consistent with dust models by Draine et al. (2007, 2014); Zubko et al. (2004).

The model is integrated with the response functions for the IRAS 100 $\mu$m and Planck 350 $\mu$m, 550 $\mu$m, and 850 $\mu$m bands, $R_{E}(\lambda)$. Since the IRAS and Planck missions assume a $\nu_{0}$/$\nu$ spectral energy distribution to compute the band flux, the models are computed using the following (with $\nu$ = $c/\lambda$):

[TABLE]

where $\nu_{\mathrm{eff}}$ is the effective frequency of the band. The corresponding wavelengths are 95.30 $\mu$m, 347.8 $\mu$m, 538.1 $\mu$m, and 830.4 $\mu$m for the IRAS 100 $\mu$m and Planck 350 $\mu$m, 550 $\mu$m, and 850 $\mu$m bands, respectively.

We pre-computed a grid of model SEDs in the IRAS and Planck bands ($S_{\mathrm{band}}^{\mathrm{mod}}$) with dust parameters $\Sigma_{d}$ $=$ $10^{-6}$ — 1 $M_{\odot}$ pc*-2* in steps of 0.02 dex, $T_{d}$ $=$ 5—49 K in steps of 0.5 K, and $\beta$ $=$ 0—3 in steps of 0.1. For each gas surface density bin, we fit the model grid to the observed, stacked SED ($S_{\mathrm{band}}^{\mathrm{mean}}$) using a probabilistic approach as in the DustBFF dust model code (Gordon et al., 2014). The dust fitting algorithm uses the base 10 logarithm of the dust surface density, since its PDF has been shown to be approximately gaussian in the log (Gordon et al., 2014), while the PDF of the dust surface density exhibits substantial wings, which could bias the estimate of the best-fit dust surface density. In other words, we compute the PDF and expectation value of $\log_{10}(\Sigma_{d}$). We compute the probability of each model as:

[TABLE]

where

[TABLE]

[TABLE]

$n$ is the number of bands (4 in our case). $\mathbb{C}$ is the covariance matrix and is the sum of the covariance matrices corresponding to calibration, background, and measurement errors, respectively: $\mathbb{C}$ $=$ $\mathbb{C}^{\mathrm{cal}}$ $+$ $\mathbb{C}^{\mathrm{bkg}}$ $+$ $\mathbb{C}^{\mathrm{meas}}$. the covariance matrix is calculated for each SED in each bin of gas surface density. The covariance matrix corresponding to measurement errors is diagonal (no covariance between bands) and the diagonal terms are given by the error on the mean dust emission in each band for a given bin of surface density: $\mathbb{C}^{\mathrm{meas}}_{ii}$ $=$ $S_{i}^{\mathrm{err}}$. The covariance matrix corresponding to residual uncertainties on the background reflects variations in the CIB and residual emission from the MW cirrus emission (the cirrus subtraction is not perfect). $\mathbb{C}^{\mathrm{bkg}}$ is estimated in regions with no detectable dust emission from the LMC and SMC. We make the same assumption as for the MW cirrus estimation in Section 3.2, and use regions where the H i emission from the LMC and SMC corresponds to H i surface densities below 4 $M_{\odot}$ pc*-2*. $\mathbb{C}^{\mathrm{bkg}}$ is given by:

[TABLE]

where $S_{i}^{k}$ is the observed surface brightness in band $i$ and pixel $k$, $N_{r}$ is the number of pixels used to estimate the background fluctuations in regions where there is no detectable dust emission ($\Sigma$(H i) $<$ 4 $M_{\odot}$ pc*-2*), $<S_{i}>$ is the mean observed surface brightness in band $i$ in those $N_{r}$ pixels, and $N_{\mathrm{stack}}$ is the number of pixels included in the stacking analysis for a given bin of gas surface density. Thus, the background covariance matrix only depends on the gas surface density bin through $N_{\mathrm{stack}}$. The other terms are constant and estimated directly in the images. In the LMC, with a quadratic fit to the cirrus emission and subtraction, $\mathbb{C}^{\mathrm{bkg}}$ is given by:

[TABLE]

In the SMC,

[TABLE]

The covariance matrix corresponding to calibration uncertainties is given by $\mathbb{C}^{\mathrm{cal}}_{ij}$ $=$ $S_{i}^{\mathrm{mean}}S_{j}^{\mathrm{mean}}$ $\bf{A}_{ij}$ and $\bf{A}_{ij}$ assumes 10%, 6.4%, 6.1%, and 0.78% absolute calibration uncertainties for the IRAS 100 $\mu$m, PLANCK 350 $\mu$m, and PLANCK 550 $\mu$m bands respectively (diagonal terms, see Planck Collaboration et al., 2016b). For IRAS 100 $\mu$m band, there is also a 2% repeatability term added in quadrature with the absolute calibration uncertainty. The correlated calibration uncertainties (non-diagonal terms) are 2% between PLANCK 350 $\mu$m and PLANCK 550 $\mu$m and 0 for the other bands. Thus,

[TABLE]

From the probabilities of each model, we compute the marginalized (1D) probability distribution functions (PDFs) for each dust parameters by integrating $P(\Sigma_{d},T_{d},\beta|S_{\mathrm{band}}^{\mathrm{mean}})$ over the other two parameters, for instance:

[TABLE]

The best dust parameters $\log_{10}\Sigma_{d},T_{d},\beta$ correspond to the expectation value of the corresponding 1D PDF. The error the parameter estimation corresponds to the square root of the variance of the 1D PDF. For instance, for the dust surface density:

[TABLE]

[TABLE]

Equations are similar for the other dust parameters. Random draws from the 1D PDFs, the best-fit model, and the data are shown in Figure 9 for the $\Sigma_{g}$ $=$ 45.9 $M_{\odot}$ pc*-2*) bin. We note that the parameter values that minimize the $\chi^{2}$ and the expectation values of each 1D PDF are consistent within errors. The SEDs and best-fit models are also shown in Figure 9 for all the gas surface density bins in the LMC and SMC.

The dust SED fitting is performed for gas surface density bins in which the mean dust emission is greater than the standard deviation in dust emission in this bin, or in other words, $S_{\mathrm{band}}^{\mathrm{mean}}$ $>$ $S_{\mathrm{band}}^{\mathrm{std}}$. The lowest gas surface density bin in which the dust can be modeled corresponds to $\Sigma_{g}$ $\sim$ 10 $M_{\odot}$ pc*-2* in both the LMC and SMC (see Figure 8). Additionally, we fit the dust SED for the dust emission obtained from a quadratic and linear estimation and subtraction of the Milky Way foreground cirrus emission.

The fractional residuals from the fits in the different bands are shown in Figure 10, and lie within 10%. This indicates that the dust SED in the different gas surface density bins is well modeled by a modified black body with a spatially varying (but spectrally constant) spectral emissivity index (the SMBB model).

3.5. Self-calibration of the dust emissivity

The dust surface density and emissivity are degenerate, since the dust emission is proportional to their product. It is therefore important to calibrate the dust FIR emissivity using our fitting technique based on independent measurements of the dust surface density. Following Gordon et al. (2014), we self-calibrate $\kappa_{0}$ by applying our dust model to the Milky Way diffuse dust SED from Compiègne et al. (2011), for which the hydrogen column density is known and the dust-to-gas ratio can be estimated from depletions derived from UV absorption spectroscopy (Jenkins, 2009). Knowing the H/D (from depletions) and gas column density (from 21 cm emission) for this sight-line, we can derive the dust surface density independently from FIR emission. The dust emissivity is then estimated by dividing the FIR emission at 160 $\mu$m by the dust surface density. We choose to calibrate the emissivity at 160 $\mu$m 1) to minimize the potential temperature dependence sometimes seen in laboratory measurements of the opacity of dust analogs at longer wavelengths (Coupeaud et al., 2011), and 2) to minimize the potential contamination by hot or stochastically heated dust. The resulting emissivity can finally be applied to the modeling of the LMC and SMC dust SED.

The diffuse MW SED in Compiègne et al. (2011), which corresponds to an H i column density $N($H i$)$ $=$ 1020 cm*-2*, has values 0.70 MJy sr*-1*, 0.55 MJy sr*-1*, 0.20 MJy sr*-1*, and 0.052 MJy sr*-1* for the 100 $\mu$m DIRBE band and 350 $\mu$m, 550 $\mu$m, 850 $\mu$m Planck bands. We have converted these values by removing the 0.77 correction for ionized gas, since the spectroscopic measurement of depletions does not include the ionized gas contribution. The corresponding uncertainties are 0.082 MJy sr*-1*, 0.048 MJy sr*-1*, 0.033 MJy sr*-1*, and 0.011 MJy sr*-1* for $N($H i$)$ $=$ 1020 cm*-2*. The H i column density toward this Galactic sight-line is $<$ 5.5$\times 10^{20}$ cm*-2*, corresponding to a median depletion factor $F_{*}$ $=$ 0.4 (see also Gordon et al., 2014). In the Milky Way, a depletion factor $F_{*}$ $=$ 0.4 corresponds to a hydrogen-to-dust ratio of $(H/D)_{\mathrm{cal}}$ $=$ 150, as obtained from the depletion patterns in Jenkins (2009).

We run this SED normalized to N(H) $=$ 1020 cm*-2* (corresponding to a hydrogen surface density $\Sigma_{\mathrm{cal}}$(H i) $=$ 0.8 $M_{\odot}$ pc*-2*) through our dust fitting code with a dummy dust emissivity of 1 cm*-2* g*-1* at 160 $\mu$m and derive a dummy dust surface density of 6.6$\times 10^{-2}$ $M_{\odot}$ pc*-2*. Since we know that the dust surface density for this SED must be $\Sigma_{d}^{\mathrm{cal}}$ $=$ $\Sigma_{\mathrm{cal}}$(H i)/$(H/D)_{\mathrm{cal}}$ $=$ 0.8/150 $=$ 5.3$\times 10^{-3}$ $M_{\odot}$ pc*-2*, the dust emissivity must be $\kappa_{0}$ $=$ 12.4 cm2 g*-1* at $\lambda_{0}$ $=$ 160 $\mu$m. The error on $\kappa_{0}$ is 50% and is dominated by errors on the dust surface density from the SED fitting. This value is consistent with different physical dust models, such as Weingartner & Draine (2001) ($\kappa_{160}$ $=$ 9.97 cm2 g*-1*), Zubko et al. (2004) ($\kappa_{160}$ $=$ 10.75—15 cm2 g*-1* depending on the fraction of graphite and amorphous carbon), Draine et al. (2007, 2014) ($\kappa_{160}$ $=$ 12.5 cm2 g*-1*), but is lower than the value of 30.2 cm2 g*-1* derived in Gordon et al. (2014) and Gordon et al. (2017). The difference might be explained by the difference in the modeling approach and the bands used. For instance, we include measurement errors, while Gordon et al. (2014) do not. We also are using four IRAS and PLANCK between 100 $\mu$m and 850 $\mu$m, while Gordon et al. (2014) use five Herschel bands between 100 $\mu$m and 500 $\mu$m. Using our fitting code with the same inputs as Gordon et al. (2014), i.e, the same Herschel 5-band SED, the same calibration covariance matrix, excluding measurement errors, we find the exact same emissivity at 160 $\mu$m (30.2 cm2 g*-1*), which validates the robustness of our fitting algorithm.

Applying the emissivity calibrated in the MW to the LMC and SMC observations assumes that the dust emissivity is similar between the MW, LMC, and SMC, which may not be correct. Indeed, Meixner et al. (2010) and Galliano et al. (2011) find evidence that the FIR emissivity of dust differs between these three galaxies. However, most of the difference in these studies lies in the spectral emissivity index, which we leave as a free parameter. Nevertheless, we also calibrate $\kappa_{160}$ based on the LMC and SMC FIR SEDs and depletion measurements corresponding to $\Sigma_{g}$ $=$ 22 $M_{\odot}$ pc*-2* in the LMC and $\Sigma_{g}$ $=$ 65 $M_{\odot}$ pc*-2* in the SMC. These gas surface densities are chosen to be in the atomic regime, since the HI-H2 transitions occur at $\Sigma_{g}$ $=$ 45 $M_{\odot}$ pc*-2* in the LMC and $\Sigma_{g}$ $=$ 110 $M_{\odot}$ pc*-2* in the SMC (Roman-Duval et al., 2014), and to minimize the measurement errors. The SEDs corresponding to those gas surface densities, shown in Figure 8, are listed in Table 2 for the LMC and SMC. Tchernyshyov et al. (2015) compiled depletion measurements in the LMC and SMC. The H/D estimates in Tchernyshyov et al. (2015) (their Figure 10) are calculated as a function of $F_{*}$, a parameter that describes the collective behavior of depletions for all elements, and is roughly correlated with mean density. To translate the variations of H/D vs $F_{*}$ into variations of G/D vs gas surface density, we 1) multiplied the H/D in Tchernyshyov et al. (2015) by 1.36 to account for Helium, and 2) examined the correlation between $F_{*}$ and $\Sigma_{g}$ and found that they are linearly correlated (with significant scatter). In the LMC F∗—$\Sigma_{g}$ relation, $F_{*}$ $=$ 0 roughly corresponds to $\Sigma_{g}$ $\sim$ 7 $M_{\odot}$ pc*-2*, while $F_{*}$ $=$ 1 corresponds to $\Sigma_{g}$ $\sim$ 43 $M_{\odot}$ pc*-2*. In the SMC, $F_{*}$ $=$ 0 corresponds to $\Sigma_{g}$ $\sim$ 6 $M_{\odot}$ pc*-2*, while $F_{*}$ $=$ 1 corresponds to $\Sigma_{g}$ $\sim$ 100 $M_{\odot}$ pc*-2*. Hence, assuming a linear relation between $\Sigma_{g}$ and $F_{*}$, the surface density of $\Sigma_{g}$ $=$ 22 $M_{\odot}$ pc*-2* chosen to calibrate $\kappa_{160}$ in the LMC corresponds to $F_{*}$ $=$ 0.42 and $\Sigma_{g}$ $=$ 65 $M_{\odot}$ pc*-2* in the SMC corresponds to $F_{*}$ $=$ 0.80. The corresponding depletion-based G/D estimates are $(G/D)_{\mathrm{dep}}$ $=$ 328 in the LMC and $(G/D)_{\mathrm{dep}}$ $=$ 800 in the SMC. Combining the estimated G/D, $\Sigma_{g}$, and FIR SEDs, we derive $\kappa_{160}$ $=$ 6.1$\pm$0.5 cm2 g*-1* in the LMC and $\kappa_{160}$ $=$ 1.65$\pm$0.2 cm2 g*-1* in the SMC. These values are lower than any estimates from dust models based on MW measurements and are consistent with the hypothesis that dust properties in the LMC and SMC differ substantially from the MW. Due to the discrepancy with physical dust models, we assume $\kappa_{0}$ $=$ 12.4 cm2 g*-1* at $\lambda_{0}$ $=$ 160 $\mu$m (based on the calibration with the MW SED) in the following, but we explore the results based on the LMC/SMC calibration in Section 4.3.

4. Results

Before getting into a quantitative analysis of the variations of the dust abundance and properties, it is worth pointing out that Figure 2 provides a clear visualization of the spatial variations of the dust abundance and properties. Bright star-forming regions appear dominantly blue/green, indicating a high ratio of bright FIR emission to H i surface density. The outskirts of the LMC and SMC and regions of low surface density (i.e., holes in the galaxies) appear to be predominantly red, indicating that the ratio of H i surface density to FIR emission is high. The color variations in Figure 2 can be explained by variations in the dust temperature (i.e., heating rate) and dust abundance (G/D). In this section, we quantify those variations using the stacking analysis and dust modeling described in Section 3.

We derive the dust surface density, dust temperature, and spectral emissivity index from the dust modeling for each interval of gas surface density following the method outlined in Section 3. The variations of these parameters are described in the following sections. Integrating the dust and gas surface densities yield total dust and gas masses summarizes in Table 3 for both MW cirrus estimation and subtraction. For comparison, Gordon et al. (2014) find total dust masses of (7.3$\pm$1.7)$\times 10^{5}$ in the LMC and (8.3$\pm$2.1)$\times 10^{4}$ in the SMC for gas masses of 4$\times 10^{8}$ for both galaxies (adjusted from the 3$\times 10^{8}$ $M_{\odot}$ H i mass to account for Helium), corresponding to G/D values of 560 for the LMC and 4900 for the SMC. Hence, we find slightly lower dust masses and higher G/D, but consistent within the errors.

4.1. Dust parameter variations with gas surface density

The resulting dust parameters are shown as a function of gas surface density in Figure 11. In both the LMC and SMC, the dust surface density increases non-linearly with gas surface density, indicating that the G/D is not constant with surface density. For comparison, a linear relation given by the slope (H/D) and intercept ($\Sigma_{I}$) quoted in Roman-Duval et al. (2014) (H/D $=$ 380${}_{-130}^{+250}$, $\Sigma_{I}$ $=$ 9.8 $M_{\odot}$ pc*-2* in the LMC and H/D $=$ 1200${}_{-420}^{+1600}$, $\Sigma_{I}$ $=$ 39 $M_{\odot}$ pc*-2* in the SMC) is shown. The H/D in Roman-Duval et al. (2014) is derived in the atomic ISM by fitting a linear function to the relation between gas and dust for $\Sigma_{H}$ $=$ 15—25 $M_{\odot}$ pc*-2* in the LMC and $\Sigma_{H}$ $=$ 40—80 $M_{\odot}$ pc*-2* in the SMC, and does not include the contribution of Helium. Thus, we multiply the H/D in Roman-Duval et al. (2014) by 1.36 to compare to the trends of G/D derived here, which include Helium in the gas surface density. Overall, the relation between dust and gas surface densities derived here is in good agreement with the linear fit derived in Roman-Duval et al. (2014) over the range where the analysis overlap. As noted before, this study allows us to derive the dust content of the LMC and SMC down to substantially lower surface densities, as shown in the top panels of Figure 11. Given that we perform a different cirrus estimation and subtraction and utilize different bands and a different fitting code, this agreement is a good test of the robustness of the results.

The SMC contains $\sim$10 times less dust than the LMC for a given gas surface density. A linear scaling with metallicity between the LMC ($Z=0.5Z_{\odot}$ Russell & Dopita, 1992) and SMC ($Z=0.2Z_{\odot}$ Russell & Dopita, 1992) would predict a factor of 2.5 difference in G/D between the LMC and SMC. This suggests that the G/D varies non-linearly with metallicity, and is consistent with the metallicity-dependence of a large sample of dwarf galaxies presented in Rémy-Ruyer et al. (2014). Chemical evolution models such as Asano et al. (2013); Feldmann (2015) have proposed explanations for this evolution: above a certain critical metallicity ($Z_{\mathrm{crit}}$ $=$ 0.1—0.2 $Z_{\odot}$ for typical star formation timescales of 0.5-5 Gyr), dust formation in the ISM is efficient enough to balance the effects of dilution by galactic inflows and outflows. Below this critical metallicity, the dust input is dominated by the non-efficient formation process in evolved stars (AGB, SN), resulting in a low dust-to-metal ratio and subsequently a low dust abundance. These models predict that, for a given star-formation timescale in the range of 0.5—5 Gyr consistent with literature (Bigiel et al., 2008; Bolatto et al., 2011; Jameson et al., 2016), the G/D in the LMC should be a factor 10 lower than in the SMC, consistent with our results.

The dust temperature occupies a small range in the LMC and SMC, with values between 25 and 30 K. There is no apparent correlation between dust temperature and gas surface density on 75 pc scales (corresponding to 5’ resolution for our dataset). Detailed 3D radiative transfer modeling would be necessary to interpret this result.

The spectral emissivity index $\beta$ appears constant within the errors across the gas surface density range, with values within 0.8—1 in the SMC and 1.2—1.4 in the LMC. Such low values of $\beta$ indicate an amorphous carbon dust grain composition (Mennella et al., 1998; Boudet et al., 2005), a conclusion also reached by Galliano et al. (2011) and Chastenet et al. (2016). In both galaxies, there is no significant correlation between $\beta$ and surface density, indicating a relatively constant ratio of amorphous carbons to silicates.

4.2. Dust abundance vs gas surface density

The gas-to-dust ratio (G/D) is computed as the ratio of the gas surface density to the dust surface density $\Sigma_{g}$/$\Sigma_{d}$. The resulting G/D is shown as a function of gas surface density in Figure 12, for the case of a linear and quadratic Milky Way cirrus emission estimation and subtraction. Errors are computed by propagating the error on $\Sigma_{d}$ described in Section 3.4 to the G/D measurement via a Monte Carlo simulation. The errors shown in Figure 12 correspond to the 50% confidence interval. We note that errors on the SED fits are dominated by calibration errors quoted in Planck Collaboration et al. (2016b), which, although not very large (10% for IRAS, 6% for Planck), can lead to substantial uncertainties on the dust temperature and therefore surface density due to the non-linear nature of the Planck function. Given the smoothness of the trends compared to the size of the error bars, it is likely that the errors are over-estimated.

We separate the gas surface density bins in two categories, ”atomic” and ”molecular”, which corresponds to regions where molecular gas fraction as traced by CO emission is lower and greater than 10%, respectively. In the SMC at 5’ resolution, molecular gas clumps are completely diluted and do not appear in the figure (i.e., over 90 pc scales, the molecular gas fraction is always $<$ 10%). This is expected due to the low metallicity of this galaxy and subsequent lack of dust shielding, leading to molecular gas filling factor that is substantially smaller than in the MW and LMC. In Figure 12, we also show the G/D vs gas surface density for the maximum possible value of $\alpha_{\mathrm{CO}}$ from Roman-Duval et al. (2014), or $\alpha_{\mathrm{CO}}$ $=$ 13 $M_{\odot}$ pc*-2* (K km s${}^{-1})^{-1}$ (3$\times$ Galactic, LMC) and 85 $M_{\odot}$ pc*-2* (K km s${}^{-1})^{-1}$ (20$\times$ Galactic, SMC). The choice of CO-to-H2 conversion factor does not affect the trends observed here, understandably because the molecular gas represents such a small fraction of the gas mass.

For comparison, the G/D derived in Roman-Duval et al. (2014) is shown in Figure 12. The G/D values derived here (black and red points) and in Roman-Duval et al. (2014) (gray band) are consistent within errors in the limited gas surface density range probed by the measurements in Roman-Duval et al. (2014). The G/D presented here is slightly higher than the value in Roman-Duval et al. (2014) for both galaxies. However, we will see in Section 5.1 that this discrepancy is due to the difference in measurement technique.Roman-Duval et al. (2014) derived the G/D as the slope (derivative) of the dust-gas relation, assumed to be linear, while we simply compute the ratio of the gas to dust surface densities in Figure 12.

Over the dynamic range observed in the LMC and SMC ($\Sigma_{g}$ $\sim$ 10—100 $M_{\odot}$ pc*-2*, shown by pink contours in Figure 5), the dust abundance (G/D) varies by a factor $\sim$ 3, from 1500 to 500 in the LMC, and by a factor $\sim$7 in the SMC from 1.5$\times$104 down to 2000. The trends are robust against the method implemented to estimate and subtract the foreground cirrus MW dust emission. Moreover, the atomic-to-molecular gas transition occurs at $\Sigma_{g}$ $=$ 45 $M_{\odot}$ pc*-2* in the LMC and $\Sigma_{g}$ $=$ 110 $M_{\odot}$ pc*-2* in the SMC (Roman-Duval et al., 2014). Therefore, most of the G/D variations observed here occur in the atomic ISM. The origins and implications of these trends are discussed in Section 5.

We note that the gas surface density, $\Sigma_{g}$, is present in the numerator of both the abscissa and ordinate in Figure 12. However, Figure 12 shows an anti-correlation between G/D and gas surface density. If the trend were due to a degeneracy between these two dependent variables, we would expect a correlation. This implies that the observed decrease of G/D with surface density is not a result of noise.

4.3. Comparison to elemental depletions

The G/D calculated from depletions (Tchernyshyov et al., 2015) is shown in Figure 12 (dashed green line). The depletion measurements are converted from H/D to G/D as explained in Section 3.5.

The slope of the G/D dependence on surface density in the depletion measurements is consistent with the slope derived here. However, the G/D derived from FIR emission and from depletions are offset by a factor of $\sim$2 in the LMC, and $\sim$8 in the SMC. A factor of 2 could easily be explained by differences in the dust properties, particularly in the FIR dust emissivity, between the LMC and the MW, where $\kappa_{160}$ is calibrated. However, a factor of 8 in the SMC is more difficult to explain with only emissivity variations. As shown in Section 3.5, reconciling the G/D measured from FIR emission with the G/D estimated from depletions would require $\kappa_{160}$ $=$ 6.1 cm2 g*-1* in the LMC and $\kappa_{160}$ $=$ 1.65 cm2 g*-1* in the SMC. This is shown in Figure 13, which shows the G/D variations in the LMC and SMC assuming those values of $\kappa_{160}$. In this case, the G/D variations with gas surface density from FIR emission and depletions agree fairly well.

While the LMC opacity value is in reasonable agreement with dust models (e.g., Draine et al., 2007, 2014; Zubko et al., 2004), the opacity calibrated on the SMC SED and depletion measurement seems unreasonably low compared to models in the literature. However, a plausible explanation comes from the comparison of the H i column density derived from UV absorption spectroscopy and from 21 cm emission in the LMC and SMC. Welty et al. (2012) find that H i column densities in the SMC derived from 21 cm emission are a factor 2—10 higher than column densities derived from Lyman $\alpha$ damping wing fitting in the surface density range probed by depletions ($\Sigma$(H i) $<$ 40 $M_{\odot}$ pc*-2*). They explain this discrepancy by the complex spatial and kinematic structure of the SMC, leading to a substantial amount of H i being located behind the stars inducing the absorption. In the LMC, the discrepancy is lower than a factor 2 owing to the much thinner structure of the disk. Thus, the differences in H i surface density between UV absorption and 21 cm measurements are consistent with the discrepancy in G/D between depletion and FIR/21 cm emission measurements observed here.

5. Discussion

5.1. G/D as a ratio or a slope

In Roman-Duval et al. (2014), we emphasized the difference between computing G/D as a ratio (integrated measurement) and as a the slope of the relation between dust and gas surface densities (local measurement). In Figure 12, there is an apparent disagreement between the G/D measurement in Roman-Duval et al. (2014) (gray band) and the study presented here (black and red points), with the G/D derived here being a factor of a few higher than in Roman-Duval et al. (2014). However, the G/D in Figure 12 is computed as a direct ratio $\Sigma_{g}$/$\Sigma_{d}$, while the Roman-Duval et al. (2014) G/D measurement shown in gray is based on a linear fit (with a non-zero intercept) to the relation between dust and gas surface densities. As suggested by the top panels of Figure 11, showing the relation between dust and gas surface densities from this analysis, the dust content of the LMC and SMC derived here is in very good agreement with the dust content observed in Roman-Duval et al. (2014). However, since the gas and dust surface densities have a non-linear (in fact super-linear) relation, their ratio is significantly different from their derivative.

In order to derive the G/D as the derivative of the dust-gas relation as in Roman-Duval et al. (2014), we fit the relation between dust and gas in the top panels of Figure 11 with different functional forms. For the LMC, the log gas-log dust relation can be approximated by a linear function (Figure 11). Therefore, a power-law (a linear fit in log-space) provides an excellent fit, which is shown as a red dashed line in Figure 11. The resulting fit in the LMC is:

[TABLE]

In the SMC, the log gas-log dust relation exhibit some curvature, which is not captured by a power-law fit. In this case, we find that a polynomial of order 3 provides the best fit to the relation between dust and gas (also shown as a red dashed line in Figure 11). The resulting relation is:

[TABLE]

We take the derivative of these functional form to derive the G/D: $G/D=d\Sigma_{g}/d\Sigma_{d}^{\mathrm{fit}}$. We bootstrap the errors on the fit using a simple Monte-Carlo simulation, by propagating the uncertainties on the $\Sigma_{d}$ measurement at each gas surface density bins into the fit and derivative. The resulting trends are shown in Figure 14 and are in excellent agreement with the G/D measured in Roman-Duval et al. (2014) in the overlapping surface density range.

5.2. Why does the dust abundance vary with surface density?

The seeds of dust grains are formed in the cool and dense envelopes of evolved stars (Bladh & Höfner, 2012) and Supernova remnants (Matsuura et al., 2011). Recent modeling and observations have shown that present-day dust abundances can only be explained if dust grows in the ISM from accretion of gas-phase metals (Boyer et al., 2012). The most direct evidence of dust growth via accretion in the ISM comes from depletions (Jenkins, 2009; Tchernyshyov et al., 2015), which indicate that the G/D decrease by a factor 3 between the most diffuse and densest sight-lines probed by UV spectroscopy of massive stars in the Milky Way, LMC, and SMC. Thus, dust grains grow in the dense ISM, and are then re-distributed to the more diffuse ISM when dense clouds are disrupted. In the diffuse ISM, dust grains are destroyed by interstellar shocks traveling at speeds 100 km s*-1* (Jones et al., 1994). Dust grains are also consumed by star-formation. In the following, we develop a simple analytical model to explain the observed trends by dust growth in the ISM.

Variations of dust abundance with density are a result of the equilibrium between production, mixing and destruction, and therefore of the differences between the corresponding timescales. Multi-phase (resolved) numerical simulations of dust evolution also find that the G/D can vary by factors of several between the diffuse and dense ISM (Bekki, 2013; Zhukovska et al., 2016; McKinnon et al., 2016). The timescale for dust formation and growth in the ISM corresponds to the collisional timescale between already existing dust grains and gas-phase metals, which in turn is inversely proportional to the gas density. Specifically, the dust growth timescale in the ISM, $\tau_{acc}$ is given in Asano et al. (2013):

[TABLE]

where $<a>$ is the mean dust grain radius (typically 0.1 $\mu$m), and the dust temperature in the LMC and SMC is $\sim$30K (see Figure 11). The metallicities of the LMC and SMC are $Z$ $=$ 0.5 $Z_{\odot}$ and $Z$ $=$ 0.2 $Z_{\odot}$, respectively (Russell & Dopita, 1992). The timescale for dust destruction by shocks is about 0.1 Gyr (Zhukovska et al., 2016, and reference therein), while star formation only consumes dust grains on $>$2 Gyr timescales (Bigiel et al., 2008). Therefore, interstellar shocks are the main dust destruction mechanism.

Observationally, elemental depletion measurements show that the dust-to-metal (D/M) ratio depends on metallicity and density. At the lowest densities ($F_{*}$ $\sim$ 0), about $f_{0}$ $=$ 20% of metals are in the dust phase in the MW and LMC, while a lower fraction (5-10%) of metals are locked in dust grains in the SMC. This number results from the timescales for dust growth, phase mixing, and dust destruction. Zhukovska et al. (2008) solve the time differential equation governing the fraction of metals in dust as a function of time:

[TABLE]

This equation also provides the density dependence of the dust-to-metal ratio (D/M) for a given time (through $\tau_{acc}$). Since abundance (M/G) are known in the LMC and SMC, we can derive the dust abundance (G/D). Therefore, we use this simple equation to model the dust abundance as a function of surface density. First, we assume that the mean density scales linearly with surface density, and that a mean density of $n_{0}$ $=$ 10 cm*-3* corresponds to a surface density of $\Sigma_{0}$ $=$ 25 $M_{\odot}$ pc*-2*:

[TABLE]

This implies that the medium has a mean density of 10 cm*-3* over a depth 100 pc, comparable to the gas disk thickness in the LMC (Elmegreen et al., 2001). Based on this relation between density and surface density, we compute $f_{d}$ as a function of surface density $\Sigma_{g}$ from Equations 15 and 16 assuming $t=10^{7}$ yr, which is roughly a dynamical time for turbulence and the timescale for dissipation of molecular clouds due to star formation. This time should be a good estimate for the timescale on which phases mix. We assume $f_{0}$ $=$ 0.2 in the LMC and $f_{0}$ $=$ 0.05 in the SMC (Tchernyshyov et al., 2015). Finally, we compute the trend in (G/D) vs $\Sigma_{g}$ via:

[TABLE]

We take $(G/D)_{0}$ to be the gas-to-dust ratio in the most diffuse regions probed in our data set, corresponding to the maximum (G/D) in the trends shown in Figure 12.

The resulting trends of (G/D) as a function of $\Sigma_{g}$ are compared to the measurements in Figure 12 (purple dashed line). Overall, we find a good agreement between this very simple model and the trends observed in the data. This is a good indication that dust indeed growth in the ISM by accretion of gas-phase metals, with density-dependent timescales consistent with simple accretion processes.

5.3. Possible effects of optically thick H i, CO-dark H2, emissivity variations, and resolution

Roman-Duval et al. (2014) pointed out the degeneracy between the effects of G/D variations, molecular gas (H2) dark in CO, optically thick H i, and variations of the dust emissivity. Dust emissivity variations and dark gas (optically thick H i and CO-dark H2) can mimic this effects of an increasing dust abundance with increasing surface density.

When dust grains coagulate in the dense ISM, their emissivity (amount of emission per unit mass) increases. Since such variations cannot be accounted for in the modeling/fitting of the dust SED due to the absence of model or observational constraints for this effect, assuming a constant emissivity implies that the surface density of dust is over-estimated and the G/D under-estimated in the dense ISM. Since coagulation only occurs at high densities, the bias only occurs in the dense ISM, leading to an apparent decrease of G/D with increasing surface density. However, this effect only kicks in at gas surface densities where a large fraction of the gas is molecular. The H i-H2 transition occurs at $\Sigma_{g}$ $=$ 45 $M_{\odot}$ pc*-2* in the LMC and 110 $M_{\odot}$ pc*-2* in the SMC (Roman-Duval et al., 2014), which is at the very high end of the surface densities probed here. The G/D variations observed here occur mostly below this threshold, and on scales larger than the typical scales of giant molecular clouds (GMCs). Hence, we do not expect coagulation and dust emissivity variations to be the cause of the trends observed here.

Optically thick H i and CO-dark H2 can lead to an apparent decrease of G/D with increasing surface density because this gas is essentially invisible to us. Roman-Duval et al. (2014) demonstrated that the atomic-to-molecular transition occurs at $\Sigma_{g}$ $\sim$ 45 $M_{\odot}$ pc*-2* in the LMC and $\Sigma_{g}$ $\sim$ 110 $M_{\odot}$ pc*-2* in the SMC. The trends observed here occur at surface densities lower than the H i-H2 transition and CO-dark H2 is therefore unlikely to affect the trends.

It is not clear how much optically thick H i contributes to the observed decrease of the G/D. However, this effect by itself could not explain the factor 4—10 variation in the G/D, or the required H i mass would be unreasonably high. Indeed, a study based on absorption 21 cm spectra in the SMC by Dickey et al. (2000b) demonstrated that the filling factor of the cold neutral phase is small ($<$15% of the H i mass) in the SMC. They derived a correction of self-absorption, which takes the form $f_{c}=1+0.667\log(N($H i$)-21.4)$ for $N($H i$)>21.4$ and $f_{c}$ $=$ 1 otherwise. At the highest column densities probed in this study ($\Sigma_{g}$ $=$ 100 corresponding to $N($H i$)$ $=$ $10^{22}$ cm*-2*), this correction factor is 1.4. Therefore, optically thick H i in the SMC does not come close to explain the decrease in G/D seen in Figure 12. No constraints on the mass of optically thick H i emission is available in the LMC, but it is unlikely that it would explain a factor 3 decrease in the observed G/D given the magnitude of the effect in the SMC.

Dust emission depends non-linearly on the dust temperature. As a result, temperature mixing in the beam can result in an underestimation of the true dust surface density is the resolution is too coarse to resolve the cold ISM (Galliano et al., 2011). However, Galliano et al. (2011) have shown that, in the LMC and SMC, the dust surface density estimate starts to suffer from this effect (at the $>$ 10% level) for resolutions coarser than 100 pc (see their Figure 6). The resolution of our maps is 75 pc (LMC) and 90 pc (SMC), and the dust surface density estimate should therefore be close to what would be derived at higher resolution. Indeed, our derived dust surface densities and mass values agree well with Gordon et al. (2014) and Roman-Duval et al. (2014), who estimated the dust content of the LMC and SMC at 10 pc resolution (Roman-Duval et al. (2014) convolved the dust surface density to 15 pc (LMC) and 45 pc (SMC), but the dust SED fitting was performed at the native Herschel resolution of 36”). The relatively low resolution of the PLANCK and IRAS data therefore does not play a role in the observed G/D trends.

5.4. Implications

The variations of the dust abundance observed in this analysis have important implications in several astrophysical fields. First, dust surface density maps or dust masses are frequently used as gas tracers, via the assumption of a constant G/D, both at low redshift (e.g., Israel, 1997; Leroy et al., 2007, 2009; Bolatto et al., 2011; Jameson et al., 2016) and high redshift (Rowlands et al., 2012). Non-linear variations of G/D with metallicity and large variations of the G/D with surface density could substantially bias the gas and dust mass estimates of galaxies.

Second, the results presented here provide some observational constraints needed by recent chemical evolution models in order to narrow down timescales for dust formation and destruction in different regions of galaxies (star-forming, diffuse, boundary with CGM). In the future, we will compare our results with some of these models in order to advance the field of dust evolution. Third, numerical models of the ISM (e.g., Glover et al., 2010; Glover & Mac Low, 2011; Smith et al., 2014) typically assume a single G/D, scaled linearly with metallicity, and constant dust properties. Since dust plays a key role in the radiative transfer in the ISM, which in turns determines the structure and composition of the different ISM phases, it is crucial to incorporate variations in the dust abundance and properties. Our results can serve as a prescription for these numerical models.

6. Conclusion

We investigate the variations of the G/D and dust properties as a function of gas surface density using a stacking analysis of H i 21 cm, CO, PLANCK and IRAS data in the LMC and SMC.

We account for systematics by correcting the H i 21 cm maps for stray-light and by carefully estimating and subtracting foreground cirrus emission from the Milky Way, using two methods (a linear and a quadratic fit to the relation between galactic H i and dust emission). We stack the LMC and SMC FIR emission observed in the IRAS 100, PLANCK 350, 550, and 850 $\mu$m bands in bins of gas surface density, and model the resulting stacked dust SED in each gas surface density interval using a modified black body fit performed within a probabilistic framework which accounts for background uncertainties, measurement, and calibration errors. We thus derive the dust surface density, temperature and spectral emissivity index as a function of gas surface density.

We find no significant variation in the dust temperature and spectral emissivity index with gas surface density on the scales probed with our all-sky surveys (75 pc in the LMC, 90 pc in the SMC). The dust temperature lies in the range 25—30 K in both galaxies and the spectral emissivity index is constant at 0.8—1 in the SMC and 1.2—1.4 in the LMC, indicating that dust grains may be pre-dominantly composed of amorphous carbon.

The dust surface density correlates non-linearly with gas surface density, resulting a G/D decreasing by factors of 4 (LMC) to 10 (SMC) between with the observed range of gas surface density ($\sim$ 10—100 $M_{\odot}$ pc*-2*). The observed trend is consistent with depletion measurements in the LMC. In the SMC, the slopes of the G/D vs gas surface density relation are consistent between depletions and our measurement, but the value of the measured G/D are discrepant by an order of magnitude, which could be explained by observational biases in the UV spectroscopic absorption H i measurements linked to the complex velocity structure of the SMC.

The G/D measurements in the LMC and SMC indicate a non-linear relation between dust abundance and metallicity, with the dust content of the SMC being a factor of 4 lower than what a linear scaling with metallicity would yield. This is consistent with the conclusions of recent chemical evolution models, which point out that dust growth in the ISM becomes too inefficient at sub-SMC metallicity to counteract metal loss by galactic-scale outflows, and that the dust content of such low metallicity galaxies should be dominated by dust input by evolved stars, resulting in a non-linear relation between G/D and metallicity.

We present a analytical model, simplified from numerical chemical evolution models, which explains the observed relation between dust and gas surface densities by accretion of gas-phase metals onto dust grains with a density-dependent timescale. This model reproduces the observations, indicating that dust growth and destruction in the ISM plays a major role in the evolution of galaxies.