DUST AND GAS IN THE MAGELLANIC CLOUDS FROM THE HERITAGE HERSCHEL KEY PROJECT. I. DUST PROPERTIES AND INSIGHTS INTO THE ORIGIN OF THE SUBMILLIMETER EXCESS EMISSION^*

The SMBB predicts the surface brightness assuming a dust population with single dust temperature modified by a simple emissivity law (Hildebrand 1983). The adopted emissivity law is

$$\begin{equation} \kappa _\lambda = \frac{\kappa _{\mathrm{eff},160}^\mathrm{S}}{160^{-\beta _\mathrm{eff}}} \lambda ^{-\beta _\mathrm{eff}}. \end{equation} \tag{ 8 }$$

The value of $\kappa _{\mathrm{eff},160}^\mathrm{S}$ is set by fitting of the diffuse MW SED (Section 5 and Table 2). The full set of fit parameters for the SMBB model are θ_S = (Σ_d, T_{eff, d}, β_eff). The values for the dust properties are effective values due to composition and temperature mixing along the line of sight and are not directly comparable to interstellar dust grain analogs studied in the laboratory (see Section 1).

The BEMBB predicts the surface brightness assuming a dust population with a single dust temperature modified by a broken emissivity law. The adopted emissivity law is

$$\begin{equation} \kappa _\lambda = \frac{\kappa _{\mathrm{eff},160}^\mathrm{BE}}{160^{-\beta _\mathrm{eff,1}}}E(\lambda) \end{equation} \tag{ 9 }$$

and

$$\begin{equation} E(\lambda) = \left\lbrace \begin{array}{@{}ll@{}}\lambda ^{-\beta _\mathrm{eff,1}} & \lambda < \lambda _b \\ \ms (\lambda _b^{\beta _\mathrm{eff,2}-\beta _\mathrm{eff,1}}) \lambda ^{-\beta _\mathrm{eff,2}} & \lambda \ge \lambda _b \\ \end{array} \right., \end{equation} \tag{ 10 }$$

where λ_b is the break wavelength and is limited to ⩾175 μm. This emissivity law is similar in form to that used by Li & Draine (2001) for astronomical silicates. The value of $\kappa _{\mathrm{eff},160}^\mathrm{BE}$ is set by the fitting of the diffuse MW SED (Section 5 and Table 2).

As we are particularly interested in measuring the submillimeter excess, we define the submillimeter excess as the excess emission at a particular submillimeter wavelength above or below that expected for a SMBB model with β_eff = β_{eff, 1}. Given the BEMBB model definition, the submillimeter excess at 500 μm is

$$\begin{equation} e_\mathrm{500} = \left(\frac{\lambda _b}{500} \right)^{\beta _\mathrm{eff,2} - \beta _\mathrm{eff,1}} - 1. \end{equation} \tag{ 11 }$$

Using e₅₀₀ as one of the fit parameters (instead of β_{eff, 2}), the fit parameters for the BEMBB model are θ_BE = (Σ_d, T_{eff, d}, β_{eff, 1}, λ_b, e₅₀₀). Note that the value of e₅₀₀ can be negative and this would indicate a submillimeter deficit. The values for the dust properties are effective values due to composition and temperature mixing along the line of sight and are not directly comparable to interstellar dust grain analogs studied in the laboratory (see Section 1).

The TTMBB predicts the surface brightness assuming two dust populations with distinctly different dust temperatures modified by a single, nonbroken emissivity law. The surface brightness is then

$$\begin{equation} S_{\lambda } = \kappa _\lambda \left[ \Sigma _{d1} B_\lambda (T_{\mathrm{eff},d1}) + \Sigma _{d2} B_\lambda (T_{\mathrm{eff},d2}) \right], \end{equation} \tag{ 12 }$$

where

$$\begin{equation} \kappa _\lambda = \frac{\kappa _{\mathrm{eff},160}^\mathrm{TT}}{160^{-\beta _\mathrm{eff}}} \lambda ^{-\beta _\mathrm{eff}}, \end{equation} \tag{ 13 }$$

the subscripts d1 and d2 refer to the two dust components, and T_{eff, d1} > T_{eff, d2}. The value of $\kappa _{\mathrm{eff},160}^\mathrm{TT}$ is set by fitting of the diffuse MW SED (Section 5 and Table 2).

For this model, the submillimeter excess at 500 μm is

$$\begin{equation} e_\mathrm{500} = \frac{\Sigma _{d2} B_\mathrm{500}(T_{\mathrm{eff},d2})}{\Sigma _{d1} B_\mathrm{500}(T_{\mathrm{eff},d1})}. \end{equation} \tag{ 14 }$$

Again, we use e₅₀₀ as a fit parameter and the full set of fit parameters for the TTMBB model are θ_TT = (Σ_d1, T_{eff, d1}, T_{eff, d2}, β_eff, e₅₀₀). Note that the value of e₅₀₀ for the TTMBB model cannot be negative, unlike the case for the BEMBB model. The values for the dust properties are effective values due to composition and temperature mixing along the line of sight and are not directly comparable to interstellar dust grain analogs studied in the laboratory (see Section 1).

It is often assumed in modified blackbody fitting that only β_eff values between one and two are valid. This is based on arguments that laboratory measurements of dust analogs only give β values between these limits. More precisely, laboratory measurements of carbonaceous and silicate dust analogs give β values between 0.8 and 2.5 for the Herschel wavelength range (e.g., Jager et al. 1998; Coupeaud et al. 2011). It is clear that luminosity weighted mixing of dust analogs with β values between 0.8 and 2.5 will always result in β_eff values in the same range. However, this is not necessarily the case for temperature mixing along the line of sight (Juvela & Ysard 2012). Combining the effects of composition and temperature mixing using full radiative transfer models, Ysard et al. (2012) give evidence that find that an β_eff (β_color in their terminology) between 0.8 and 2.5 is reasonable for a range of realistic cases. Thus, we include versions of the SMBB, BEMBB, and TTMBB models that have β_eff values restricted to be between 0.8 and 2.5. However, we caution that it is more statistically correct to include β_eff values outside this range as measurement noise can create SEDs that require nonphysical β_eff values to provide statistically robust fits.

Our models produce SEDs that are well sampled in wavelength, but our observations have a very coarse wavelength sampling as they are taken through filters with broad response functions. It is important to correctly model the effects of these broad response functions on the models to give accurate fits to the observations. For this paper, we start with the model predictions of the surface brightnesses at a wavelength resolution that well resolves the PACS and SPIRE bandpasses (Müller et al. 2011b; Griffin et al. 2013). Then, the band surface brightnesses were determined by integrating over their respective band response functions using

$$\begin{equation} S_\mathrm{band} = \frac{ \int S_\nu R_E(\nu) d\nu }{\int (\nu _o/\nu)^{-1} R_E(\nu) d\nu }, \end{equation} \tag{ 15 }$$

where R_E(λ) is the response function appropriate for extended sources given in fractional transmitted energy. The ν_o = c/λ_o values are given by λ_o = 100, 160, 250, 350, and 500 μm for the bands with the same names. Equation (15) mathematically models the data that is produced by the PACS and SPIRE instruments and data reduction pipelines. The integration is done in energy units (e.g., MJy sr⁻¹) as both instruments use bolometers that measure energy (not photons). The denominator of this equation normalizes R_E(λ) and accounts for the PACS and SPIRE calibration convention where the calibration is given at specific wavelengths (λ₀) and for a S(ν) = ν⁻¹ reference spectrum.

For this work, the covariance matrix is defined as

$$\begin{equation} \mathbb {C} = \mathbb {C}_\mathrm{cal} + \mathbb {C}_\mathrm{bkg}, \end{equation} \tag{ 19 }$$

where $\mathbb {C}_\mathrm{cal}$ is the absolute surface brightness covariance matrix and $\mathbb {C}_\mathrm{bkg}$ is the background covariance matrix. The units of these covariance matrices are (MJy sr⁻¹)².

The $\mathbb {C}_\mathrm{cal}$ is given by the details of the PACS and SPIRE absolute flux calibrations. The SPIRE instrument has been calibrated using a model of Neptune with an absolute uncertainty correlated between bands for point sources of 4% and a repeatability that is uncorrelated between bands of 1.5% (Griffin et al. 2013; Bendo et al. 2013). For extended sources, it is recommended to add an additional 4% to account for the correlated uncertainty in the total beam area resulting in an 8% correlated uncertainty between bands (Herschel Space Observatory 2011). The PACS instrument has been calibrated using models of stars and asteroids with an absolute uncertainty correlated between bands for point sources of 5% and a repeatability uncorrelated between bands of 2% (Müller et al. 2011a; Balog et al. 2013). Similar to SPIRE, for extended sources we add an additional 5% correlated uncertainty to account for uncertainties in the total beam area resulting in a 10% correlated uncertainty between bands. Finally, we assume the PACS and SPIRE calibrations are independent given that PACS is calibrated using stars and SPIRE using Neptune. Given this information, the elements of $\mathbb {C}_\mathrm{cal}$ are

$$\begin{equation} (\mathbb {C}_\mathrm{cal})_{ij} = S^\mathrm{mod}_i(\theta) S^\mathrm{mod}_j(\theta) \left[ (\mathbf {A}_\mathrm{cor})_{ij} + (\mathbf {A}_\mathrm{uncor})_{ij} \right], \end{equation} \tag{ 20 }$$

where

$$\begin{equation} \mathbf {A}_\mathrm{cor} = \left[\begin{array}{@{}ccccc@{}}0.1^2 & 0.1^2 & 0 & 0 & 0 \\ \ms 0.1^2 & 0.1^2 & 0 & 0 & 0 \\ \ms 0 & 0 & 0.08^2 & 0.08^2 & 0.08^2 \\ \ms 0 & 0 & 0.08^2 & 0.08^2 & 0.08^2 \\ \ms 0 & 0 & 0.08^2 & 0.08^2 & 0.08^2\end{array} \right] \end{equation} \tag{ 21 }$$

and

$$\begin{equation} \mathbf {A}_\mathrm{uncor} = \left[ \begin{array}{@{}ccccc@{}}0.02^2 & 0 & 0 & 0 & 0 \\ \ms 0 & 0.02^2 & 0 & 0 & 0 \\ \ms 0 & 0 & 0.015^2 & 0 & 0 \\ \ms 0 & 0 & 0 & 0.015^2 & 0 \\ \ms 0 & 0 & 0 & 0 & 0.015^2\end{array} \right]. \end{equation} \tag{ 22 }$$

The background covariance matrix, $\mathbb {C}_\mathrm{bkg}$ is calculated empirically from a large set of pixels visually identified as lying outside of the emitting region of each galaxy. The background pixels are in the full images and were processed as described in Section 2. The terms of the covariance matrix are calculated using

$$\begin{equation} \sigma _{ij}^2 = \frac{ \sum _k^N (S_i^k - \langle S_i \rangle) (S_j^k - \langle S_j \rangle) }{N-1}, \end{equation} \tag{ 23 }$$

where N is the number of background pixels, $S_i^k$/$S_j^k$ is the ith/jth band of the kth pixel, and 〈S_i〉/〈S_j〉 is the average background in the ith/jth band. For the LMC, N = 52113 and

$$\begin{equation} \mathbb {C}_\mathrm{bkg}({\rm LMC}) = \left[ \begin{array}{@{}ccccc@{}}4.23 & 0.78 & 0.65 & 0.33 & 0.14 \\ \ms 0.78 & 2.37 & 0.85 & 0.43 & 0.18 \\ \ms 0.65 & 0.85 & 0.91 & 0.47 & 0.20 \\ \ms 0.33 & 0.43 & 0.47 & 0.25 & 0.11 \\ \ms 0.14 & 0.18 & 0.20 & 0.11 & 0.057\end{array} \right] \end{equation} \tag{ 24 }$$

and for the SMC, N = 4012 and

$$\begin{equation} \mathbb {C}_\mathrm{bkg}({\rm SMC}) = \left[ \begin{array}{@{}ccccc@{}}2.64 & 0.56 & 0.30 & 0.14 & 0.064 \\ \ms 0.56 & 1.18 & 0.46 & 0.23 & 0.094 \\ \ms 0.30 & 0.46 & 0.36 & 0.20 & 0.089 \\ \ms 0.14 & 0.23 & 0.20 & 0.12 & 0.054 \\ \ms 0.064 & 0.094 & 0.089 & 0.054 & 0.030\end{array} \right]. \end{equation} \tag{ 25 }$$

These empirical covariance matrices illustrate that background is highly correlated with the correlation increasing in strength toward longer wavelengths. This is illustrated by the correlation matrix (terms are $\mathbb {C}_{ij}/[\sigma _i \sigma _j]$) for the SMC:

$$\begin{equation} \mathrm{corr}_\mathrm{bkg}({\rm SMC}) = \left[ \begin{array}{@{}ccccc@{}}1.00 & 0.32 & 0.31 & 0.25 & 0.23 \\ \ms 0.31 & 1.00 & 0.70 & 0.61 & 0.49 \\ \ms 0.30 & 0.70 & 1.00 & 0.94 & 0.85 \\ \ms 0.25 & 0.61 & 0.94 & 1.00 & 0.91 \\ \ms 0.23 & 0.49 & 0.85 & 0.91 & 1.00\end{array} \right]. \end{equation} \tag{ 26 }$$

The LMC correlation matrix is very similar and thus is not shown. The positive and nonzero correlation terms are signatures that the correlated noise in the background is due to real astronomical signals. In this case, it is traceable to the residual foreground MW cirrus emission and the integrated emission from background galaxies. The higher covariance values for the LMC is a reflection of the increased difficulty of background subtraction for this galaxy.

The fitting technique we use fully computes the nD likelihood function that a particular model fits the SED of a pixel where n is the number of fit parameters. One way to visualize the results is to create 1D likelihood functions for each fit parameter by marginalizing (integrating) over all the other parameters. This is shown in Figure 1 for the BEMBB model for a single pixel in the SMC for three different assumptions; assuming uncorrelated uncertainties, including the full covariance, and including the full covariance while restricting β_{eff, 1} to vary between 0.8 and 2.5. The results for pixels in the LMC are similar. With the same overall uncertainties, we obtain a much narrower function with a stronger likelihood by including the known covariance between the bands (Section 4.1) than by assuming that there is no correlation between bands. In this case, including the known covariance between bands results in better constraints on the fit parameters as the allowed model space is reduced. The impact of a limited β_{eff, 1} range is shown in this figure where, not surprisingly, it makes for a narrower 1D likelihood function than allowing β_{eff, 1} to vary in order to fully sample the β_{eff, 1} 1D likelihood function. Note that this limitation simply crops the β_{eff, 1} 1D likelihood function, but changes the shape of the other 1D likelihood functions significantly.

Zoom In Zoom Out Reset image size

The goal of the sensitivity tests is to determine if there are systematic shifts in recovered parameters and if the uncertainty on the recovered parameters matches that measured from the widths of the 1D likelihood functions. We simulated observations by picking a model SED and adding noise using the Cholesky factorization of the covariance matrix appropriate as if the model was observed like the SMC was observed in HERITAGE. The results using the LMC noise model are very similar. We repeated the simulation for each model SED 20 times to provide a good sampling of the recovered fit parameter uncertainties and systematic offset from the input fit parameters.

As we are testing the ability of this fitting technique to recover parameters by fitting simulated observations, this requires a way to measure the recovery of the input model parameters. The main output of the fitting is the nD likelihood function, but it is often useful to distill these results to the"best-fit" or summary values. We use three different ways to define the "best-fit" values. The first is the most traditional definition of the "best fit" and corresponds to the maximum likelihood ("max"). This is also called the "traditional χ²" method in some papers (e.g., Kelly et al. 2012; Juvela et al. 2013). The "max" value is most useful when plotting the best-fitting model with observations or investigating the fitting residuals. The second is the expectation value ("exp") which is the likelihood weighted average of the parameter and is a reflection of the full likelihood function. This "exp" value reflects the best "average" value as it reflects the full likelihood function (not just the peak like the "max" value). We find the "exp" particularly useful for making images of the fit parameters. The third way to reflect the best fit is take a realization of the full nD likelihood function itself ("realize"). This involves randomly sampling the likelihood function and reflects the full likelihood function's shape in a statistical sense. The "realize" method is most useful when studying the ensemble behavior of the fit parameters for many pixels.

The results for runs with 2000 randomly picked BEMBB models are shown in Figure 2. All three different methods of determining the "best-fit" parameters give similar results with similar trends with each parameter. The "exp" gives the lowest systematic error in the recovery, but the "max" gives the lowest scatter. The "realize" method provides a nominally worse recovery than both the other methods, but is a fuller picture of true sensitivity of the fitting. Overall, which "best-fit" method used depends on the particular question being asked. We illustrate this later in this paper and in the companion paper on the gas-to-dust ratio (Roman-Duval et al. 2014).

Zoom In Zoom Out Reset image size

Of particular interest for this paper is the fact that the recovery of the submillimeter excess, e₅₀₀, is good to around 10%, on average, for the "realize" method and around 1% for the "exp" method. For the companion paper (Roman-Duval et al. 2014), the fit parameter of main interest is Σ_d and the recovery is good, on average and in log (Σ_d) units, to 0.05 for the "realize" method and 0.001 for the "exp" method. This excellent recovery of log (Σ_d) holds even in the presence of significant scatter in T_{eff, d} and may be due to other parameters in the fitting varying to compensate. Note that, for the "exp" method, we computed the expectation value of log (Σ_d) because we found that the sensitivity tests showed significantly less systematic bias than if we computed the expectation value of Σ_d. We confirmed that the widths of the 1D likelihood functions match the noise in the recovery of the input model parameters.

The number of parameters in our models are three, five, and five for the SMBB, BEMBB, and TTMBB models, respectively. In this paper, these models are fit to FIR-submillimeter SEDs that are composed of five data points. At first glance, this violates the rule that fitting requires at least one data point more than the number of fitting parameters to provide a unique solution. This is correct, if the fitting is done with a model that can fit any distribution of data points. This is clearly not the case for our models since they are all constrained to have a spectral shape of one or two modified blackbodies. In other words, they cannot fit arbitrary spectral shapes, but are constrained by our knowledge of the physics of dust grain emission. Effectively, we are using more than just five data points in our fits because we combine the data points with a larger body of observations that informs our understanding of dust physics and, therefore, the appropriate models to use. Finally, our use of full likelihood functions explicitly accounts for the impact of the number of parameters on how well we can determine each fit parameter. Using full likelihood functions has the additional benefit of measuring how well each parameter is constrained by the data explicitly. Some parameters are better constrained than others as shown in Figure 1. For example, Σ_d and T_{eff, d} are better constrained as the overall level and spectral shape are well constrained by the observations, but the detailed spectral shape is less well constrained and this strongly impacts β_{eff, 1}, λ_b, and e₅₀₀.

For the diffuse MW emission, we use the Compiègne et al. (2011) measurement where emission was measured by correlating the IR versus H i emission maps in atomic gas dominated regions of the MW. The IR measurements we use are mainly the COBE/FIRAS spectrophotometry from 127 to 1200 μm supplemented by the DIRBE 100 μm photometry. Because we want to calibrate our models using the same bands as used for the HERITAGE observations, we integrated this diffuse MW SED using the method described in 3.5 for all the bands except the PACS 100 μm band. For this band, we adopted the DIRBE 100 μm measurement since the bandpasses are similar. The resulting MW diffuse SEDs are 0.71, 1.53, 1.08, 0.56, and 0.25 MJy sr⁻¹ (10²⁰ H atom)⁻¹ for the 100, 160, 250, 350, and 500 μm and are plotted in Figure 3. These values differ from those given for the same bands by Compiègne et al. (2011) mostly because we have not included the 0.77 correction for ionized gas. In addition, there are minor differences in the response curves used. We do not include the 0.77 correction for ionized gas because the depletion measurements do not include any ionized gas correction. For the uncertainties, we have assumed 5% correlated and 2.5% uncorrelated terms (see Section 3.5) given the high quality of the COBE, FIRAS, and DIRBE calibrations.

Zoom In Zoom Out Reset image size

Because the MW diffuse SED is measured as a correlation between dust and gas emission, the constraint we need is the MW diffuse gas-to-dust ratio. We use the work of Jenkins (2009) to determine the appropriate gas-to-dust ratio since this work provides an excellent compilation and summary of MW depletions. The observed H columns of our adopted FIR-submillimeter MW diffuse SED are log [N(H)] < 20.7. The average depletion of all the sightlines with these column densities tabulated by Jenkins (2009) is F_* = 0.36. F_* is the depletion factor and measures the overall depletions in a sightline. Using the depletion fits of Jenkins (2009) with F_* = 0.36, the diffuse MW gas-to-dust ratio is computed to be 150.

We calibrate the value of κ_{eff, 160} in each of our models so that they reproduce the MW diffuse observed gas-to-dust ratio of 150. For our work, we have chosen 160 μm to set our normalization of κ_{eff, λ} because shorter wavelengths have a weaker dependence on temperature based on laboratory investigations of dust analogs (Coupeaud et al. 2011). The κ_{eff, 160} values required for each model based on the "exp" method of determining the best fits (see Section 4.3) are given in Table 2. The second uncertainty on κ_{eff, 160} is an estimate of the systematic uncertainty (see next paragraph). The fit parameters for each model are also given in this table, along with 1σ uncertainties. The larger relative uncertainties on κ_{eff, 160} for the BEMBB model as compared to the SMBB can be directly traced to the larger number of BEMBB fit parameters. The "max" best-fit models are plotted in Figure 3.

The κ_{eff, 160} values for the SMBB and BEMBB models agree favorably with other determinations while the value for the TTMBB model does not. For example, if "astronomical" silicate grains with a = 0.1μm and ρ = 3 g cm⁻³ are used, then κ_{eff, 160} = 13.75 cm² g⁻¹. Such grain properties are often assumed for simple modified blackbody fits because this is the average size for a Mathis et al. (1977) grain size distribution (Hildebrand 1983). The widely used Weingartner & Draine (2001) full dust grain model for R(V) = 3.1 has a κ_{eff, 160} = 9.97 cm² g⁻¹. The updated version of this model has a κ_{eff, 160} = 12.5 cm² g⁻¹ (Draine & Li 2007; Draine et al. 2014). The κ_{eff, 160} values for the Zubko et al. (2004) models that include graphite and amorphous carbon range from 10.75 to 15.0 cm² g⁻¹. Finally, the Weingartner & Draine (2001) model for the SMC Bar extinction curve with no 2175 Å extinction feature has κ_{eff, 160} = 13.1 cm² g⁻¹. Using the range of these model κ_{eff, 160} values, we estimate that there is a ±2.5 cm² g⁻¹ additional uncertainty on κ_{eff, 160} due to systematic uncertainties in our knowledge of dust grains.

The TTMBB model with κ_{eff, 160} = 517 ± 214 cm² g⁻¹ requires a dust grain that is very efficient at emission, yet this level of efficiency is much higher than any astronomically reasonable dust grain. A much simpler explanation is that the dust in the MW diffuse ISM is not well modeled by a TTMBB model that includes a very cold (T_{eff, d} ∼ 6 K) dust grain population. This is the same conclusion given by the Reach et al. (1995) analysis of the FIRAS data. There still may be regions in the ISM of the MW or other galaxies that are well described by the TTMBB model. To allow for such regions, we adopt the κ_{eff, 160} of the SMBB model as the value for the TTMBB model.

The variations in the κ_{eff, 160} values in the literature and between the different models used in this paper clearly indicate that κ_{eff, 160} is sensitive to the model assumptions. Thus, it is important to calibrate each model explicitly with the diffuse MW SED and a depletion measured gas-to-dust ratio. This is a standard calibration method for dust grain models (Draine & Li 2007; Compiègne et al. 2011) and we advocate that such calibrations be done for all dust emission models (Bianchi 2013). Such model calibrations will allow for meaningful comparisons between the results from different models.

One obvious question is: which model, SMBB, BEMBB, or TTMBB, fits the observations best? The answer to this question will give an indication of the origin of the submillimeter excess. The most straightforward method to test how well a model fits the data is to examine the residuals of the data to the fits. The χ² value computed using Equation (18) gives such a quantitative measure of the residuals. For the SMC, the pixel averaged χ² value is 3.47 for the SMBB model, 0.88 for the BEMBB model, and 1.83 for the TTMBB. The models with 0.8 < β_eff < 2.5 have higher average χ² values than the unconstrained versions. For example, the β_eff constrained version of the BEMBB model for the SMC has an average χ² value of 1.32. The LMC average χ² values behave similarly.

More evidence that the BEMBB fits the data best (out of the three models) can be found by examining the behavior of the fit residuals versus surface brightness. Figure 4 shows the fit residuals for the SPIRE 250 μm band for all three models used in this paper for both Magellanic clouds. The trends for other bands are similar, especially in the relative behavior of the fit residuals between the models. This figure clearly shows that the simplest model (SMBB) has residuals larger than expected given the known uncertainties. This holds for β_eff unconstrained and constrained to be between 0.8 and 2.5. In addition, the residuals for the SMC have a systematic trend with more negative residuals at intermediate surface brightnesses. Such a trend is not consistent with the uncertainties in the absolute flux calibration or the background subtraction. Of all models, the BEMBB model without any constraint on β_eff fits the data best. Overall, the BEMBB model shows the smallest residuals with no obvious trend with surface brightness unlike the other models. The BEMBB model consistently shows smaller residuals in all the bands, not just the SPIRE 250 μm band. The other models have higher overall residuals and show systematic offsets and/or trends with surface brightness. The BEMBB and TTMBB models have the same number of fit parameters, yet the behavior of their residuals are different. This illustrates that it is not only the number of fit parameters that is critical for the fitting accuracy, but also the allowed spectral shapes.

Zoom In Zoom Out Reset image size

Overall, the BEMBB spectral shapes fit the data better than the TTMBB and SMBB spectral shapes. This is evidence that the submillimeter excess is more likely to be due to emissivity variations than a second population of cold dust.

The total dust masses are of interest for studies of the lifecycle of dust in the LMC and SMC (Boyer et al. 2012; Matsuura et al. 2013; Zhukovska & Henning 2013). In addition, they can be used along with the total gas masses as a way to tell if a model produces realistic amounts of dust (see Section 6.3).

We give the dust masses for the different models in Table 3 integrated over the >3σ pixels. The restricted β_eff version of the models produces results that are very similar and are not given in the table. The dust mass values are given as total ± statistical uncertainty ± uncertainty due to the κ_{eff, 160} uncertainty. To convert from dust surface density to dust mass, we use distances of 60 kpc (Hilditch et al. 2005) and 50 kpc (Walker 2012) for the SMC and LMC, respectively. The total dust masses are computed from the "realize" method to produce dust surface density maps that provide a full accounting of the likelihood functions for all pixels. Ten different maps were made for each galaxy using the "realize" method that samples the likelihood function once for each pixel. This provides a robust measurement of the impact of the fitting noise of each pixel in the integrated dust mass measurement. The average and statistical uncertainty of the integrated dust mass were computed from the 10 maps. The large number of pixels in each galaxy result in the total dust mass changing only slightly between different realizations and this is the origin of the small statistical uncertainty. These dust masses are integrated only over the areas that were detected at 3σ above the background in all five Herschel bands measured by HERITAGE. Pixels >3σ contribute 0.79, 0.73, 0.62, 0.61, and 0.61 of the SMC global fluxes of 15.7, 20.8, 14.5, 8.3, and 3.9 kJy for the PACS100, PACS160, SPIRE250, SPIRE350, and SPIRE500, respectively. For the LMC, these fractions are 0.91, 0.89, 0.87, 0.87, and 0.87 for global fluxes of 223, 259, 142, 73, and 31 kJy for the same bands. The global fluxes quoted here differ from those given by Meixner et al. (2013) due to our subtraction of MW cirrus foreground and the additional background subtraction step.

The quantitative impact of correctly including the correlated noise in the measurements can be illustrated by noting that assuming the noise is uncorrelated between bands results in the BEMBB model giving fits with a total SMC dust mass that is ∼50% higher than the total dust mass given in Table 3. The importance of accounting for the full likelihood function is equally important: the total SMC dust mass for the BEMBB model is ∼50% higher using the "max" values and ∼30% lower using the "exp" values of log (Σ_d) when compared to the "realize" value given in Table 3. The "realize" values are the correct values for determining the total dust mass values as they statistically reflect each pixel's full likelihood function, asymmetries and all, in the sum of the individual pixel masses. The "max" and "exp" values only reflect a limited portion of the likelihood function and this systematically biases the results. This is additional evidence that the likelihood functions for Σ_d are not well behaved Gaussians centered on the "max" value (see Figure 1).

Our total dust masses are only lower limits because we do not include the dust responsible for the emission with surface brightnesses below 3σ in any band. We can estimate the dust mass due to these <3σ regions by modeling the integrated flux of these regions for each galaxy. Basically, we fit the SED that is the difference from the global fluxes quoted above and the integrated fluxes from <3σ pixels. The resulting integrated dust masses for the BEMBB model and the <3σ pixels are (5.9 ± 3.6) × 10⁴ and (1.6 ± 1.3) × 10⁴ M_☉ for the LMC and SMC, respectively. The uncertainties are quite large due to the low surface brightnesses and strong mixing of environments in these integrated SEDs. Combining the <3σ pixel dust masses with those for >3σ pixel (Table 3), we find total dust masses of (7.3 ± 1.7) × 10⁵ and (8.3 ± 2.1) × 10⁴ M_☉ for the LMC and SMC, respectively. For reference, the total gas masses that correspond to the same areas and same background removal as these total dust masses are 3.1 × 10⁸ and 3.0 × 10⁸ M_☉ for the LMC and SMC, respectively.

Bot et al. (2010b) obtained global dust masses for both galaxies by fitting Draine et al. (2007) dust models to their global fluxes. They found masses of 3.6 × 10⁶ and 0.29–1.1 × 10⁶ M_☉ for the LMC and SMC, respectively. Leroy et al. (2007b) fit the spatially resolved Spitzer observations with (Dale & Helou 2002) models and find a total SMC dust mass of 3 × 10⁵ M_☉. These values of the dust masses are factors of four to five larger than our values. The differences are likely due to different assumptions in the models used, the fitting techniques, the broader wavelength range of data, and/or the increased mixing of environments.

One test of the submillimeter excess origin is to investigate how the overall gas-to-dust ratios for each model compare to the expected ratios. We explore overall gas-to-dust ratios as a test of the consistency of each dust model with expectations based on the measured gas masses and metallicities of the LMC and SMC. The detailed spatial behavior of the gas-to-dust ratio with the environment is investigated in Roman-Duval et al. (2014).

The gas-to-dust ratios for each galaxy and all three models are given in Table 3. The dust masses are integrated over all the pixels that are detected at >3σ in all observed bands. The total H gas masses given in the table footnote are integrated for the same pixels as the dust masses. The H i masses are directly from the H i measurements (Stanimirović et al. 2000; Muller et al. 2003; Kim et al. 2003) without any correction for opaque H i (Dickey et al. 2000; Fukui et al. 2014). The H₂ masses are computed from CO observations (Mizuno et al. 2001, 2006; Fukui et al. 2008; Wong et al. 2011) using X_CO = 4.7 × 10²⁰ (Hughes et al. 2010) for the LMC and X_CO = 6 × 10²¹ (Leroy et al. 2007a) for the SMC. The appropriate X_CO to use is a matter of debate, but the expected range of this conversion factor is not large enough to strongly impact the total gas masses (Fukui & Kawamura 2010; Bolatto et al. 2013). The ratios given only include hydrogen, and thus are formally H gas-to-dust ratios, but for simplicity we refer to them as gas-to-dust ratios.

The range of reasonable gas-to-dust ratios can be estimated in three ways. The first scales the range of observed gas-to-dust ratios in the MW by the LMC and SMC metallicities. The second assumes the MW depletion factors and applies them to the measured LMC and SMC abundances. The third assumes that all the metals available are in the form of dust and this produces a minimum possible gas-to-dust ratio. The MW depletions and gas-to-dust ratios vary with environment and the global values in the Magellanic clouds will be some unknown mix of different ISM environments. As a result, we can only predict a possible range of gas-to-dust ratios.

The first method assumes that the relative amount of metals in the LMC and SMC dust is the same as the MW, but scaled in proportion to each galaxy's metallicity. Thus, the expected gas-to-dust ratio will be two times (LMC) and five times (SMC) the MW gas-to-dust ratio. The MW gas-to-dust ratio varies from ∼250 for the very diffuse ISM (F_* = 0) to ∼100 for the moderately dense ISM (F_* = 1) (Jenkins 2009). For the LMC, we therefore expect a gas-to-dust ratio between 200 to 500 while, for the SMC, we expect a gas-to-dust ratio between 500 and 1250. The second method assumes the MW depletion patterns (Jenkins 2009) and the measured LMC and SMC abundances for each element (Russell & Dopita 1992). The resulting expected LMC gas-to-dust ratios range between 150 to 360 and the expected SMC gas-to-dust ratios range between 540 to 1300. Combining the two different methods, the expected gas-to-dust ratios are 150–500 and 500–1300 for the LMC and SMC, respectively. Finally, the minimum allowed gas-to-dust ratio can be computed by assuming all the metals in the ISM in the form of dust. Assuming the measured LMC and SMC abundances, this gives minimum gas-to-dust ratios of 105 and 300, respectively. These expected gas-to-dust ratios are given in Table 3.

The gas-to-dust ratios for all three models are plotted in Figure 5 along with the allowed ranges for reasonable depletions and maximum depletions. From Table 3 and this figure, it is clear that the TTMBB models give gas-to-dust ratios that are lower than even possible assuming that all the metals are present in the dust. The TTMBB model gives low gas-to-dust ratios because it requires large dust masses for the second cold component to be able to reproduce the observed submillimeter excess emission. Thus, the TTMBB model is not a reasonable model for the dust emission in the LMC or SMC. The SMBB and BEMBB models give similar gas-to-dust ratios for both galaxies. For the LMC, both models give ratios that are well within the reasonable range of values. For the SMC, these two models both give values that are above the reasonable values. This is an indication that the depletions in the SMC are lower than those the in MW or that the dust properties are different (e.g., a smaller κ_{eff, 160} value than that assumed in this paper).

Zoom In Zoom Out Reset image size

The spatial variations across both galaxies in the different fit parameters for the BEMBB model are shown in Figures 6 and 7. We only show the BEMBB results here because the evidence in the previous sections gives a fairly strong indication that the BEMBB fits the data best (Section 6.1) and provides a physically reasonable gas-to-dust ratio (Section 6.3). The maps of dust surface density (Σ_d) and temperature (T_{eff, d}) show qualitatively similar behaviors to previous works (Bot et al. 2004; Leroy et al. 2007b; Bernard et al. 2008). In detail, our maps differ mainly in showing finer structure due to the higher spatial resolution Herschel observations. One illustration of this effect is that the peak T_{eff, d} in the 30 Dor region in our map is ∼60 K, significantly higher than the ∼35 K found by Bernard et al. (2008).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The higher spatial resolution of our maps does allow for detailed investigations of individual star-forming regions. This is illustrated by Figure 8 where cutouts of the BEMBB fit parameter maps for a star-forming region in each galaxy are shown. The morphology of these two star-forming regions is similar. The SPIRE 250 μm emission is strongly peaked in the region centers in contrast to the dust surface density which is more constant across the regions. This difference is caused by the center of these regions having high T_{eff, d} values. The β_eff and e₅₀₀ maps of both regions have very similar morphologies, visually illustrating that these two fit parameters are strongly correlated. Finally, the λ_b images show coherent structures with fairly small variations overall. The submillimeter excess as parametrized by e₅₀₀ is near zero in the center of the two star-forming regions and rises rapidly to values around one near the edges. This behavior is intriguing, but the strong correlations of e₅₀₀ with β_eff indicate that more work is needed to determine if this is real or due to noise induced correlations.

Zoom In Zoom Out Reset image size

The overall properties of the global submillimeter excess between the LMC and SMC show trends that are consistent with previous work. The average LMC and SMC e₅₀₀ values are 0.27 and 0.43 when the average is done using the "realize" method and each pixel has equal weight. This can be visually seen in the e₅₀₀ images in Figures 6 and 7 where the SMC shows a higher filling factor of high e₅₀₀ values than the LMC. This trend of the lower metallicity SMC having a higher submillimeter excess than the LMC is expected given the results from global studies of the submillimeter excess Rémy-Ruyer et al. (2013). A fairer comparison of the absolute value of e₅₀₀ with global SED fits is the dust surface density weighted averages that are 0.11 and 0.26 for the LMC and SMC, respectively. Finally, the average values of λ_b are ∼240 for both types of averages and both galaxies. This wavelength is similar to that found by Li & Draine (2001) from fitting the DIRBE MW diffuse spectrum.

To investigate the variations in fit parameters more quantitatively, we plot all the correlations between the different fit parameters for the LMC in Figure 9. The plots for the SMC are very similar and are not shown. These plots show the density of points where each point represents a single pixel. The values used for each pixel use the "realize" method where the likelihood functions are randomly sampled once for each pixel. This means that these density plots statistically sample the full information for the fit from each pixel. Repeating the "realize" method process with a different random sampling for each pixel produces plots that are very similar. This indicates that these plots fully capture the correlations between fit parameters with a single sampling of each pixel's likelihood function due to the large number of pixels. Plots created using the "max" and "exp" methods are significantly different because they do not fully include the information on the uncertainties in the fits to each pixel. As an example of the difference between the different "best-fit" methods, a flat likelihood function would show a single value for "max" and "exp," while the "realize" method would have a value that was randomly distributed over the entire parameter range.

Zoom In Zoom Out Reset image size

These plots show that many of the parameters are correlated with each other, sometimes quite strongly. The strongest correlations are seen between log (Σ_d) and T_{eff, d}, log (Σ_d) and β_{eff, 1}, T_{eff, d} and β_{eff, 1}, and β_{eff, 1} and e₅₀₀. The origin of these correlations can be either real or a result of interactions between noise in the measurements and model fit parameters. The correlation between log (Σ_d) and T_{eff, d} is real in that it reflects the detection thresholds of the HERITAGE data. Hotter dust can be more easily detected at lower dust surface densities than cooler dust due to the $T_\mathrm{eff,d}^4$ behavior of blackbodies. The anti-correlation between T_{eff, d} and β_eff is one of the correlations that has been studied extensively to learn if it is due to noise or real variations in the dust properties (Dupac et al. 2003; Shetty et al. 2009a, 2009b; Galliano et al. 2011; Juvela & Ysard 2012; Kelly et al. 2012; Ysard et al. 2012; Veneziani et al. 2013). Laboratory data on dust analogs do show a shallow anti-correlation between T_{eff, d} and β_eff (Coupeaud et al. 2011), but noise in measurements also produces a similar or larger anti-correlation. Kelly et al. (2012) have proposed the use of a hierarchical Bayesian model to solve for the true T_{eff, d}–β_eff correlation, where the hierarchical model assumes a single T_{eff, d} and β_eff with some distribution around these values. In fitting an entire galaxy, such an assumption is not justified because, for example, there are regions near star formation that will be significantly hotter than regions further away. In addition, Juvela et al. (2013) find there are biases in all the currently proposed methods for determining the true T_{eff, d}–β_eff relation. Thus, we choose to graphically display the correlations using the "realize" method and not explicitly fit for the correlation. In future work, we plan to incorporate additional observations of the ISM and physical models for the correlations between different ISM parameters (e.g., dust and gas surface densities).

Figure 9 shows the correlations between the submillimeter excess e₅₀₀ and other dust properties. The value of e₅₀₀ is positively correlated with Σ_d and β_{eff, 1} and negatively correlated with T_{eff, d}. This may be real or it may be due to the T_{eff, d} versus β_{eff, 1} anti-correlation that is also very clearly seen. The positive correlation between e₅₀₀ and Σ_d is the opposite of what was found by Galliano et al. (2011) for a pathfinder study using a portion of the HERITAGE data on the LMC and Paradis et al. (2012) for the MW. The difference between these works and our work may be due to changes in the PACS and SPIRE calibration, different fitting methods, and/or different dust emission models. Future work will investigate these differences by using the same data, same fitting code, and expanding the dust emission model to include more sophisticated dust emission models.