Benchmarking Dust Emission Models in M101

We correct the images for extended-source calibration appropriately and convert them to the same units in the following way. We apply extended-source correction factors for the IRAC images. These are multiplicative corrections of 0.91, 0.94, 0.695, and 0.74 at 3.6, 4.5, 5.8, and 8.0 μm, respectively, as suggested by the IRAC Instrument Handbook. ¹⁰ For SPIRE, we adopt calibration factors of 90.646, 51.181, and 23.580 MJy/sr/(Jy/beam) at 250, 350, and 500 μm, respectively, as suggested by the SPIRE Handbook. ¹¹ These factors convert the data to MJy sr⁻¹ from Jy/beam and also correct from point-source to extended-source calibration. The PACS data in units of Jy/pixel are converted to units of MJy sr⁻¹ using the image pixel size (contained in the headers).

We then process all of the maps following the same steps, as described below. We remove a background in each image at their native resolution by fitting a 2D plane using regions identified not to have significant galaxy emission. We find these regions with the following procedure: (1) we use the r₂₅ radius (${r}_{25}\,\sim 12^{\prime}$) to first mask the galaxy (this covers the visible SPIRE 500 emission completely); (2) we measure the median of the remaining pixels and the standard deviation, and clip all of those that are three standard deviations above that median; (3) we iterate that clipping until the medians at iterations i and i + 1 differ by less than 1%. This clipping is only done for the purpose of measuring a background level, and we do not keep this mask applied to the data for the following steps. Table 2 lists the final standard deviations of the background pixels in each band.

Each map is then convolved to the SPIRE 500 PSF (FWHM ∼ 36'') using convolution kernels from Aniano et al. (2011). Finally, all of the maps are aligned and projected onto the astrometric grid of the SPIRE 500 image. In Figure 1, we show the 16 bands that we use to model the dust emission from 3.4 to 500 μm.

Zoom In Zoom Out Reset image size

The final pixel size (9'') oversamples the SPIRE 500 beam size, to which all data are convolved. We take this into account by correcting by $\sqrt{{N}_{\mathrm{pix}}/{N}_{\mathrm{beam}}}$ when calculating uncertainties on quantities that use the values of multiple pixels.

In Section 6, we use gas surface density and metallicity measurements in addition to IR SEDs. We use the same data as Chiang et al. (2018) to get the total gas surface density, Σ_gas, combining gas surface density maps from H i 21 cm and CO(2 → 1) emission converted to H₂. The H i 21 cm line data is from the THINGS collaboration (The H i Nearby Galaxy Survey; Walter et al. 2008). The map is converted from integrated intensity to surface density assuming optically thin 21 cm emission, following Walter et al. (2008, Equations 1 and 5, and multiplying by the atomic mass of hydrogen).

We use CO (2 → 1) emission from HERACLES (Leroy et al. 2009) and convert it to H₂ surface density, assuming a line ratio R₂₁ = 0.7 and a α_CO conversion factor, in units of ${M}_{\odot }/{\mathrm{pc}}^{2}\,{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}$, following two prescriptions (both including helium): the first one is representative of the MW CO-to-H₂ conversion, ¹²

$$\begin{eqnarray}&&{\alpha }_{\mathrm{CO}}^{\mathrm{MW}}=4.35,\end{eqnarray} \tag{ 1 }$$

and the second one follows Bolatto et al. (2013),

$$\begin{eqnarray}&&\,\begin{array}{l}{\alpha }_{\mathrm{CO}}^{\mathrm{BWL}13}=2.9\ {e}^{(0.4/Z^{\prime} )}{\left({{\rm{\Sigma }}}_{\mathrm{Total}}^{100}\right)}^{-\gamma }\\ \mathrm{with}\left\{\begin{array}{cc}\gamma =0.5\ & \mathrm{if}\ {{\rm{\Sigma }}}_{\mathrm{Total}}^{100}\geqslant 1\\ \gamma =0\ & \mathrm{otherwise}\end{array}\right.,\end{array}\,\end{eqnarray} \tag{ 2 }$$

where $Z^{\prime}$ is the metallicity relative to solar metallicity, ¹³ ${{\rm{\Sigma }}}_{\mathrm{Total}}^{100}$ is the total surface density map in units of 100 M_⊙ pc⁻². The gas maps are convolved to a 41'' Gaussian PSF (Chiang et al. 2021); because we plot the radial profile of the gas maps by averaging in growing annuli, we find that the minor differences between the dust and gas maps resolutions are negligible. We measure an uncertainty of ∼ 0.3 Kkm s⁻¹ for the CO(2 → 1) map and an uncertainty of ∼0.4 M_⊙ pc⁻² for the H i 21 cm map. ¹⁴ We build a radial profile by averaging Σ_gas in growing annuli from the center out to r₂₅.

Metallicity measurements $12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$ are taken from the CHAOS survey (Berg et al. 2020), derived from 72 H ii regions of M101. In particular, we use their fitted metallicity gradient

$$\begin{eqnarray}&&12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})=8.78\pm 0.04-0.75\pm 0.07\ (r/{r}_{25})\end{eqnarray} \tag{ 3 }$$

to convert it to a radial profile of metal surface density (see Section 6.1). We adjusted the slope given by Berg et al. (2020, 0.90) to account for the different r₂₅ used here. Berg et al. (2020) find that the scatter in individual region metallicities around the measured gradient is ∼10%.

Tracing the total metal mass from nebular emission lines of oxygen is subject to systematic uncertainties. As pointed out by Jenkins (2009), oxygen shows unexplained depletion patterns, compared to how much is expected to be locked in the solid phase. However, Peimbert & Peimbert (2010) measured the depletions of heavy elements in Galactic and extragalactic H ii regions and found minimal depletion of oxygen (≲0.1 dex). Other elements may be of use for tracing metallicity (as is the case in Berg et al. 2020), but given the widespread use of O/H for extragalactic studies, we use the typical $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ tracer.

We use single-temperature modified blackbody models. For these two models, we only fit photometry from 100 to 500 μm. At λ < 100 μm, stochastically heated grains contribute to the emission and are not well modeled by a modified blackbody. Assuming the optically thin case, the dust emission I_ν is described as

$$\begin{eqnarray}&&{I}_{\nu }(\lambda )={\kappa }_{\nu }(\lambda )\ {{\rm{\Sigma }}}_{\mathrm{dust}}\ {B}_{\nu }(\lambda ,{T}_{\mathrm{dust}}),\end{eqnarray} \tag{ 4 }$$

where B_ν (λ, T_dust) is the Planck function at wavelength λ (in MJy sr⁻¹), T_dust the dust temperature, Σ_dust the dust surface density, and κ_ν the opacity. We use the simple power-law opacity (Equation (5)) and a broken power-law opacity (Equation (6)), which we normalize at λ = 160 μm. The broken power-law model presented the best results in terms of quality of fits in the study by Chiang et al. (2018) and yielded physically reasonable Σ_dust and T_dust values in their study.

We follow Gordon et al. (2014) and Chiang et al. (2018) and calibrate the opacity values for each model. We fit the modified blackbody model to the dust emission per H column of the MW high-latitude cirrus described in Chiang et al. (2018) to derive the opacity. The abundance constraints are based on a depletion strength factor typical for lines of sight with N_H similar to that of the MW cirrus, e.g., F_* = 0.36 (Jenkins 2009). This sets the allowed dust mass per H atom. By fitting the temperature for the MW cirrus, we can then derive the opacity for each modified blackbody model. More details on the opacity and comparison to the physical models can be found in Section 3.2.

3.1.1. Simple Emissivity

In this case, the opacity is described as a single power law:

$$\begin{eqnarray}&&{\kappa }_{\nu }={\kappa }_{{\nu }_{0}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{\beta },\end{eqnarray} \tag{ 5 }$$

where β is the spectral index. For all blackbody models, we fix λ₀ = 160 μm and ν₀ = c/λ₀. The free parameters for this model are then the dust surface density Σ_dust, the dust temperature T_dust, and the dust spectral index β. In this model, the calibrated ${\kappa }_{{\nu }_{0}}=10.10\pm 1.42\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1}$ (Chiang et al. 2018) from fitting the high-latitude cirrus as described above.

3.1.2. Broken Emissivity

The broken-emissivity model stemmed from the identification of the submillimeter excess (Gordon et al. 2014). It allows a change of the dust spectral index with wavelength, meant to better reproduce the far-IR slope than does a simple modified blackbody. In this model, the value of the opacity is wavelength dependent, such that

$$\begin{eqnarray}{\kappa }_{\nu }=\left\{\begin{array}{ll}{\kappa }_{{\nu }_{0}}{\left(\displaystyle \frac{\nu }{{\nu }_{0}}\right)}^{\beta } & \lambda \lt {\lambda }_{c}\\ {\kappa }_{{\nu }_{0}}{\left(\displaystyle \frac{{\nu }_{c}}{{\nu }_{0}}\right)}^{\beta }{\left(\displaystyle \frac{\nu }{{\nu }_{c}}\right)}^{{\beta }_{2}} & \lambda \geqslant {\lambda }_{{\rm{c}}}\end{array}\right.,\end{eqnarray} \tag{ 6 }$$

where λ_c is the wavelength at which the opacity changes and ν_c its equivalent frequency. Following Chiang et al. (2018), we fix β = 2 and λ_c = 300 μm. The free parameters are then the dust surface density Σ_dust, the dust temperature T_dust, and the second dust spectral index β₂. The calibrated opacity, ${\kappa }_{{\nu }_{0}}$, for this model is 20.73 ± 0.97 cm² g⁻¹ (Chiang et al. 2018), as described above.

The physical dust models used in our study have opacities that are set by each individual model calibration procedure. All models use similar constraints from the MW high-latitude cirrus (described in more detail in Appendix A). For ease of comparison to the opacities we have derived for the modified blackbody models (10.1 cm² g⁻¹ for the simple-emissivity and 20.7 cm² g⁻¹ for the broken-emissivity model at 160 μm), we list the opacity in the physical models below all scaled to 160 μm for comparison:

• 1.  
  Draine & Li (2007): from the Weingartner & Draine (2001b) model updated in DL07, we report κ₁₆₀ = 10.2 cm² g⁻¹ (see also Bianchi 2013);
• 2.  
  Compiègne et al. (2011): from the work of Bianchi (2013), we report κ₁₆₀ = 12.0 cm² g⁻¹;
• 3.  
  Jones et al. (2017, THEMIS): from the work of Galliano et al. (2018), we report κ₁₆₀ = 14.2 cm² g⁻¹;
• 4.  
  HD21: the models use κ₁₆₀ ∼ 9.95 cm² g⁻¹ (B. Hensley 2021, private communication).

For all models, we use the listed opacity regardless of the specific environment in M101 we are studying. Few of the currently available physical dust models have been calibrated in any environment other than the MW cirrus (although the Draine & Li 2007 model does have Small and Large Magellanic Cloud–like calibrations, though they are not widely used even for these very galaxies; Sandstrom et al. 2010; Chastenet et al. 2019). Indeed, it is standard in current extragalactic applications to apply MW cirrus R_V = 3.1 dust calibrations across all environments (e.g., Davies et al. 2017; Aniano et al. 2020). In detail, this is unlikely to be correct because the opacity can and probably should evolve as a function of the environment, and it is clear that a single F_* = 0.36 value does not describe the depletion in the ISM over the full range of column densities probed in galaxies. However, for the purposes of our comparative study of widely used dust models, we proceed by using the R_V = 3.1 MW cirrus calibrations from each model. Even if there were a potential way to adjust opacity with the H column in M101, work by Ysard et al. (2015) using THEMIS suggests that this is insufficient to predict changes in dust properties.

Physical dust models assume a composition, density, and shape for the dust grains and adopt heat capacities and optical properties from laboratory and theoretical studies that are appropriate for such materials. For simplicity, most models assume spherical grains or planar molecules for PAHs. In the case of PAHs, the grains/molecules are additionally described by an ionization state, which changes the absorption cross-sections as a function of wavelength. The temperature and emission of a grain of a given size, shape, and composition in a radiation field with a specified intensity and spectrum can then be calculated analytically (e.g., Draine & Lee 1984; Desert et al. 1986).

The full dust population is represented by a grain size distribution and abundance relative to H for the specified compositions. Physical dust models are calibrated by adjusting the grain size distributions and dust mass per H to simultaneously match observations of extinction, emission, and abundances (and more recently, polarization) in a location where the underlying radiation field that is heating the dust is well known. This has generally been taken to be the high-latitude MW cirrus, where the radiation field intensity and spectrum are approximately given by the Mathis et al. (1983) model for the solar neighborhood. The degree to which the models must adhere to the somewhat uncertain abundance constraints varies from model to model. For example, the modified blackbody models from Gordon et al. (2014) use the depletion measurements in the MW as a strict limit. However, most physical models allow the final element abundances to vary from depletions, using them only as a loose guide.

In the following, we use four physical dust models: Draine & Li (2007), Compiègne et al. (2011), THEMIS (Jones et al. 2017), and HD21. Here we briefly describe these models, the key differences between them. Details on their respective calibration methodologies can be found in Appendix A.

3.3.1. Draine & Li (2007)

3.3.2. Compiègne et al. (2011)

3.3.3. THEMIS

3.3.4. B. S. Hensley & B. T. Draine (2021, in preparation)

Because each model was developed independently, it is not always possible to create model grids that have exactly the same parameter sampling, because of the lack of parameter equivalence. However, we attempt to do so as much as possible.

Radiation field—We implement identical dust heating in each of the physical models: a fraction γ of the dust mass in a pixel is heated by a power-law distribution of radiation fields with ${U}_{\min }\lt U\leqslant {U}_{\max }$ (Dale et al. 2001), where U is the interstellar radiation field, expressed in units of the MW diffuse radiation field at 10 kpc from Mathis et al. (1983). The remaining fraction of dust (1 − γ) is heated by the minimum radiation field ${U}_{\min }$ (Draine & Li 2007):

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{1}{{M}_{\mathrm{dust}}}\displaystyle \frac{{{dM}}_{\mathrm{dust}}}{{dU}} & = & (1-\gamma )\ {U}_{\min }\\ & & +\gamma \ \displaystyle \frac{(\alpha -1)}{{U}_{\min }^{1-\alpha }-{{\rm{U}}}_{\max }^{1-\alpha }}\ {U}^{-\alpha }.\end{array}\,\end{eqnarray} \tag{ 7 }$$

We fix α = 2 for all models; ${U}_{\max }={10}^{7}$ for the DL07, THEMIS, and MC11 models; and ${U}_{\max }={10}^{6}$ for HD21. This choice is constrained by the available parameter ranges in the models: the ${U}_{\max }$ values for the DL07 and HD21 models are fixed, and we do not have the freedom to change them; thanks to DustEM , we can adjust ${U}_{\max }$ for MC11 and THEMIS. ¹⁹ For THEMIS and the MC11 model, which have multiple dust populations that can be independently heated, each population is heated by the same radiation field.

The minimum radiation field parameter grid is fixed by the values provided in the Draine & Li (2007) model. Thanks to the DustEM tool, we can use the exact same values with THEMIS and the MC11 model. These parameter values in HD21 however are not exactly identical to those in ${U}_{\min }^{\mathrm{DL}07}$. We therefore use the spectra with the closest ${U}_{\min }^{\mathrm{HD}21}$ values. Although not strictly equal, the radiation field values in Hensley & Draine are within 5% of ${U}_{\min }^{\mathrm{DL}07}$.

For the modified blackbody models, we use a single radiation field intensity and translate it into a dust temperature. We use the relationship

$$\begin{eqnarray}&&U\propto {T}_{\mathrm{dust}}^{\beta +4}\quad \mathrm{with}\ \beta =2\end{eqnarray} \tag{ 8 }$$

to convert ${U}_{\min }$ to T_dust and find the approximately matching sampling to use in the blackbody models. We use the normalization from Draine et al. (2014), i.e., U = 1 at T_dust = 18 K, found using the same radiation field from Mathis et al. (1983).

Mid-IR feature carriers—The q_PAH parameter is kept strictly identical between DL07 and HD21. We choose to use THEMIS and the MC11 models in a similar fashion. Despite the definition of "PAHs" being different in THEMIS, we parameterize the model so that the fraction of small grains can vary, keeping the large-carbonaceous-to-silicate grains ratio constant. The MC11 default model has four populations that can vary independently, and we choose to tie together SamC, LamC, and aSil to leave the amount of PAHs as a free parameter. We explain these choices in more detail in Section 7.

Stellar surface brightness—In addition to the dust model parameters, we use a scaling parameter of a 5000 K blackbody, Ω_*, to model the stellar surface brightness visible at the shortest wavelengths. This temperature is a good approximation as the shortest wavelengths are nearly on the Rayleigh–Jeans tail. The free parameter Ω_* scales the amplitude of the stellar blackbody.

In Table 3, we list the final free parameters we use for each model. There are five free parameters for the physical models, and three for the modified blackbody models. Figure 2 shows all of the dust emission models used in this study. The top-left panel shows the fiducial MW high-galactic-latitude diffuse ISM models, labeled "Galactic SED." The IR MW diffuse emission from Compiègne et al. (2011) is also plotted. The bottom-left panel shows the physical dust models at the same radiation field ${U}_{\min }$ = 1. The other different panels detail the models: THEMIS is divided into two grain populations when fitted to the SED of M101 (note that we tie the lCM20 and aSilM5 population); for the MC11 model, we tie the SamC, LamC, and aSil emissions together, and thus there are two free parameters for the dust mass: f_PAH or f_sCM20, and the total dust surface density. The two modified blackbodies are shown with their best-fit values to the MW diffuse SED: {T_dust = 20.9 K, β = 1.44} for the simple-emissivity model, and {T_dust = 18.0 K, β₂ = 1.55} for the broken-emissivity model.

Zoom In Zoom Out Reset image size

4.2.1. Bayesian Fitting with DustBFF

We use the DustBFF tool from Gordon et al. (2014) to fit the data with the chosen dust models. DustBFF is a Bayesian fitting tool that uses flat priors for all parameters. The probability that a model with parameters θ fits the measurements (S^obs) is given by

$$\begin{eqnarray}&&{ \mathcal P }({{\boldsymbol{S}}}^{\mathrm{obs}}| \theta )=\displaystyle \frac{1}{Q}{e}^{-{\chi }^{2}(\theta )/2},\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}\begin{array}{rcl}{Q}^{2} & = & {\left(2\pi \right)}^{{\rm{n}}}\ \det | {\mathbb{C}}| \\ {\chi }^{2}(\theta ) & = & {\left[{{\boldsymbol{S}}}^{\mathrm{obs}}-{{\boldsymbol{S}}}^{\mathrm{mod}}(\theta )\right]}^{{\rm{T}}}\ {{\mathbb{C}}}^{-1}\ [{{\boldsymbol{S}}}^{\mathrm{obs}}-{{\boldsymbol{S}}}^{\mathrm{mod}}(\theta )],\end{array}\end{eqnarray} \tag{ 10 }$$

where S ^X is the observed (X = obs) or modeled ($X=\mathrm{mod}$) 16 band SEDs used here, and ${\mathbb{C}}$ is the covariance matrix that includes uncertainties from random noise, astronomical backgrounds, and instrument calibration (described further below).

To create ${{\boldsymbol{S}}}^{\mathrm{mod}}(\theta )$, each model spectrum is convolved with the spectral responses for the photometric bands used here. PACS and SPIRE band integrations are done in energy units, whereas all the others are done in photon units, as necessitated by the instrument's calibration scheme. We also follow the reference spectra used for the calibration of each instrument: the MIPS bands require a reference shape in the form of the 10⁴ K blackbody while the other bands use a reference shape as 1/ν.

4.2.2. Covariance Matrices

The covariance matrix in the previous equations describes the uncertainties on the measured flux, both due to astronomical backgrounds, and instrumental noise and the uncertainties in calibrating the instruments. This takes into account the correlation between the photometric bands due to the calibration scheme and the correlated nature of astronomical background signals. We define a background matrix and an instrument matrix such that ${\mathbb{C}}={{\mathbb{C}}}_{\mathrm{bkg}}+{{\mathbb{C}}}_{\mathrm{ins}}$ to propagate the correlated errors of the background and noise, and of the calibration uncertainties, respectively.

Background covariance matrix ${{\mathbb{C}}}_{\mathrm{bkg}}$—The "background" in our images encompasses astronomical signals that do not come from the IR emission of the target M101. It is dominated by different objects depending on the wavelength and can therefore be correlated between bands. For instance, the background from 3.4 to 5.8 μm is mostly from foreground stars, as well as zodiacal light; from 8 to 24 μm, it is from evolved stars and background galaxies; in the far-IR, it is dominated by Galactic cirrus emission and background galaxies. To include this uncertainty, we measure this combination of signals in "background pixels" using the processed data and a masking procedure described in Section 4.2.3.

The elements of the background covariance matrix are calculated as

$$\begin{eqnarray}&&{\left({{\mathbb{C}}}_{\mathrm{bkg}}\right)}_{i,j}^{2}=\displaystyle \frac{{\displaystyle \sum _{k}}^{{\rm{N}}}({S}_{i}^{k}-\langle {S}_{i}\rangle )({S}_{j}^{k}-\langle {S}_{j}\rangle )}{{\rm{N}}-1},\end{eqnarray} \tag{ 11 }$$

where S_X ^k is the flux of pixel k in band X, and 〈S_X 〉 is the average background emission in band X (close to 0). The final number of 9'' pixels used to measure the covariance matrix is N = 18,920.

We display the background correlation matrix in Figure 3 (not the elements' absolute values but the Pearson correlation coefficients; see Gordon et al. 2014). Three clear sets of positive correlations appear, as previously explained: correlations are due to starlight in the near-IR, evolved stars and MW cirrus in the mid-IR, and background galaxies and MW cirrus in the far-IR. Additionally, some bands may be noise dominated: it is the case for the PACS 70 band, which is only weakly correlated with other bands.

Zoom In Zoom Out Reset image size

Instrument calibration matrix ${{\mathbb{C}}}_{\mathrm{ins}}$—This matrix is calculated as the quadratic sum of the correlated and uncorrelated errors for each instrument. The correlated errors (matrix with diagonal and anti-diagonal terms) refer to the instrument calibration itself (m_cal in Table 2, while the uncorrelated errors (matrix with diagonal terms only) express the instrument stability, or repeatability (m_sta in Table 2). We use the same calibration errors reported in Chastenet et al. (2017, see Table 2) for the Spitzer /MIPS and Herschel bands. The correlated errors for the IRAC bands were changed to take into account the uncertainties in the extended-source correction factors, larger than the calibration error themselves. The errors due to repeatability are unchanged. The errors for WISE bands were taken from the WISE documentation. ²⁰

The elements of the calibration matrix ${{\mathbb{C}}}_{\mathrm{ins}}$ are calculated "model by model" as

$$\begin{eqnarray}&&{\left({m}_{\mathrm{ins}}\right)}_{i,j}^{2}={S}_{i}^{\mathrm{mod}}(\theta )\ {S}_{j}^{\mathrm{mod}}(\theta )\ ({m}_{\mathrm{cal},i,j}^{2}+{m}_{\mathrm{sta},i,j}^{2}),\end{eqnarray} \tag{ 12 }$$

with particular elements

$$\begin{eqnarray*}\begin{array}{rcl}{m}_{\mathrm{cal},i,j} & = & 0\ \mathrm{if}\ i,\,j\ \mathrm{belong}\ \mathrm{to}\ \mathrm{two}\ \mathrm{different}\ \mathrm{instruments};\\ {m}_{\mathrm{sta},i,j} & = & 0\ \mathrm{if}\ i\ne j.\end{array}\end{eqnarray*}$$

We fit all the pixels that are above 1σ of the background values, in all bands. We use these pixels to show parameter maps and radial profiles, while the galaxy-integrated values are calculated for pixels above 3σ detection above the background (black contours in Appendices B and C).

4.2.3. Stars in the Background/Foreground

Figure 4 shows an example of the fits for each model in a single pixel (marked by the cross in Figure 1). In each panel, we plot the best-fit spectrum (colored lines) to the data SED (empty circles). The residuals are shown by the colored symbols. Negative residuals mean that the model overestimates the data. For example, the negative (and decreasing) residuals from 250 to 500 μm in the physical model panels are representative of a systematic overestimation of the data (present also in other locations; see Appendix C, Figures 19–22). In Figure 5, we show the fractional residuals (Data-Model)/Model in each band for the pixels above the 3σ threshold. Appendix D shows the reduced χ² for all models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bulk of the residuals in the short-wavelength bands are within the instrument uncertainties and calibration errors. For example, despite the larger uncertainty due to the extended-source correction, the IRAC 8 band shows residuals mostly within 10%. Below 4 μm, THEMIS shows a clear offset that may be related to the absence, or low amount, of an ionized component in the HAC population, which leads to enhancement of the mid-IR features. It is also worth noting that these bands are dominated by starlight, modeled by a 5000 K blackbody, which is independent of the dust model itself. The residuals at 12 μm are systematically positively offset by less than 10%, but all physical models show very narrow residual distributions. This is in contrast with the broader residuals in the IRAC 8 and WISE 22 bands, where we can see more differences between models.

All models show differences in the central values of the residual distribution at 100 μm. At 160 μm, all residuals overlap, and models perform fits of similar quality. At longer wavelengths, significant differences begin to appear. At all far-IR wavelengths, the modified blackbody models reproduce the SEDs the best. This is likely because of the additional parameter that can adjust the far-IR slope (β in the simple-emissivity model and β₂ in the broken-emissivity model). The SPIRE 250 band shows symmetrical residuals centered on 0 for the physical models (except THEMIS), while the residuals get progressively worse at 350 and 500 μm for all physical models, showing that on average, the modeled far-IR slope of the SEDs is steeper than the data.

In the SPIRE 350 and SPIRE 500 bands, the large number of pixels underestimated by the models show residuals much larger than the uncertainties, ruling out statistical noise and indicating that the models are not able to fit these wavelengths. In the SPIRE 500 band, some of the pixels show the so-called "submillimeter excess" seen in other studies (e.g., Galametz et al. 2014; Gordon et al. 2014; Paradis et al. 2019).

We compute several galaxy-averaged quantities: the total dust mass M_dust, the dust mass-weighted average radiation field $\langle \overline{U}\rangle$ for the physical models, and the mass-weighted average dust temperature 〈T_dust〉 for the blackbody models. The average radiation field $\overline{U}$ is calculated for each pixel as

$$\begin{eqnarray}&&\overline{U}=(1-\gamma )\ {U}_{\min }+\gamma \times {U}_{\min }\ \displaystyle \frac{\mathrm{ln}({{\rm{U}}}_{\max }/{U}_{\min })}{1-{U}_{\min }/{U}_{\max }},\end{eqnarray} \tag{ 13 }$$

because we fixed α = 2 (Aniano et al. 2020). The galaxy-integrated mass-weighted averages are calculated as

$$\begin{eqnarray}\,\begin{array}{rcl}\langle \overline{U}\rangle & = & \displaystyle \frac{{\displaystyle \sum }_{j}{\overline{U}}_{j}\times {{\rm{\Sigma }}}_{\mathrm{dust},j}}{{\displaystyle \sum }_{j}{{\rm{\Sigma }}}_{\mathrm{dust},j}}\\ \langle {T}_{\mathrm{dust}}\rangle & = & \displaystyle \frac{{\displaystyle \sum }_{j}{T}_{\mathrm{dust},j}\times {{\rm{\Sigma }}}_{\mathrm{dust},j}}{{\displaystyle \sum }_{j}{{\rm{\Sigma }}}_{\mathrm{dust},j}}.\end{array}\,\end{eqnarray} \tag{ 14 }$$

The integrated values are calculated over the pixels above the 3σ detection threshold. In Figure 6, we show these measurements for each model as diagonal elements. We also provide the ratios between all models in M_dust and $\langle \overline{U}\rangle$ (or 〈T_dust〉) to explicitly show their differences. These off-diagonal elements read as Y-model/X-model (e.g., ${M}_{\mathrm{dust}}^{\mathrm{HD}21}/{M}_{\mathrm{dust}}^{\mathrm{DL}07}=0.77$).

Zoom In Zoom Out Reset image size

The top panel of Figure 6 shows the total dust masses. The broken-emissivity modified blackbody model yields the lowest total dust mass, while the MC11 model yields the highest. The simple-emissivity model shows a very different spatial variation of dust surface density from the other models (see Section 5.3). In a recent study of the KINGFISH sample, Aniano et al. (2020) found a total dust mass of 9.14 × 10⁷ M_⊙ (before their renormalization) by fitting the DL07 at MIPS 160 resolution (∼39''). ²³ Using the CIGALE SED fitting tool (e.g., Boquien et al. 2019 and references therein) and THEMIS, Nersesian et al. (2019) found a total dust mass of 4.70 × 10⁷ M_⊙ , which is similar to the 5.05 × 10⁷ M_⊙ from THEMIS, in this study. It is worth noting that Nersesian et al. (2019) performed fits to the integrated SED, which could lead to a lower dust mass (Aniano et al. 2012; Utomo et al. 2019). However, low signal-to-noise pixels are included in the integrated SED and excluded in the resolved fits.

It is interesting to note that despite the fairly good agreement (∼10%, Figure 5) of all physical models with each other (and modified blackbody models at far-IR wavelengths) in reproducing the data, the differences in dust masses can be much larger. This suggests intrinsic opacity values rather than fit quality dominate the differences between models in dust mass. This is supported by the recent extensive study done by Fanciullo et al. (2020). By comparing the literature opacities (including three models used in this work) with laboratory dust analog opacities, they found that dust masses can be overestimated by more than an order of magnitude.

The bottom panel of Figure 6 shows the integrated values for $\langle \overline{U}\rangle$ and 〈T_dust〉. As expected, the HD21 model requires the highest radiation field, based on its colder dust. THEMIS and the MC11 model show similar values of $\langle \overline{U}\rangle$ . Nersesian et al. (2019) found a dust temperature of 21.7 K by fitting a modified blackbody (with the THEMIS opacity), close to that yielded by the modified blackbody models here.

The mass-weighted average temperatures correspond to radiation fields of 2.5 and 2.3 (using Equation (8)) for the broken-emissivity and the simple-emissivity models, respectively. They are in good agreement with the radiation field values fitted by the physical models.

Figure 7 shows maps of the dust surface density, Σ_dust, for each model, as well as their ratios with each other. We can see that the physical dust models DL07, HD21, THEMIS, and MC11 are all fairly close to each other (light colors), with variations in dust surface density within a factor of 2. They all yield similar dust surface density structures and appear to vary from one another by a spatially smooth offset. THEMIS and the HD21 model show the closest Σ_dust values but show an inversion of their ratio around 0.38 r/r₂₅, where THEMIS requires less dust (see Figure 12). The HD21/DL07 and MC11/THEMIS ratio maps are particularly flat, with ratios of ∼1.3–1.4 in both cases, pointing at the resemblance in their large grains properties and size distribution. The MC11 model requires the most dust mass. This is consistent with the comparison analysis in Chastenet et al. (2017). The dust surface density from the broken-emissivity model is consistently lower than the physical dust models. It shows a rather smooth offset, which indicates that the spatial variations are fairly similar between them. The dust surface density from the simple-emissivity model shows different spatial structures in the center and the outskirts of the galaxy. It yields high Σ_dust values in the center but drops more rapidly than any other model with increasing radius (Figure 12; see also Chiang et al. 2018).

Zoom In Zoom Out Reset image size

A caveat of our approach is the difference in treating the heating of dust grains between physical dust models and modified blackbodies. In the latter, we only use a single temperature, which is not equivalent to the ensemble of radiation fields used in the physical models. However, given the two extreme behaviors of the modified blackbodies (the simple-emissivity model requiring a high dust mass, and the broken-emissivity model requiring the lowest), it is not obvious that the use of multiple temperatures (or radiation fields) drives the differences in dust mass observed here. Rather, the calibration of the modified blackbodies, and their effective opacity, seems to be more important.

We perform the same ratio analysis with $\overline{U}$ , derived from the fitted parameters in the physical dust models (Equations 13 and 14). In Figure 8, we show the radial profiles for $\overline{U}$ , as well as the parameter maps and the ratios of each model. In the radial profile, the thick lines stop where the selection effect due to fitting only bright pixels becomes important. Using the SPIRE 500 image, we found that the radial profile of IR emission for all pixels and for fitted pixels only differ significantly at ∼6' (i.e., 0.5 r₂₅). The variations in $\overline{U}$ are reflective of those of ${U}_{\min }$, because the γ values are overall small, lending more power to the delta function than the power law (Equation (7)).

Zoom In Zoom Out Reset image size

The overall variations of $\overline{U}$ appear to be rather smooth, which is expected because it is dominated by the diffuse radiation field ${U}_{\min }$. However, we do find enhanced values of $\overline{U}$ in H ii regions. The ratio maps do not strongly show these peaks, which indicates that all models behave similarly and require higher radiation field intensities in these regions. As for the Σ_dust parameter, the ratio maps for $\overline{U}$ do not display any conspicuous spatial differences, but rather an offset between each model. This is also visible in the upper panel of Figure 8.

The HD21 model shows the highest values of $\overline{U}$. This is expected as the dust in the model is "colder" than other models, due to very few large carbonaceous grains. This leads to a higher radiation field intensity required to reach the same luminosity.

The spatial distributions of the γ parameter are similar in all physical models, and we chose not to show them. Instead, we use the average radiation field, $\overline{U}$, which includes γ in its calculation. The HD21 model shows the lowest values of γ, which, combined with the highest values of ${U}_{\min }$, means it requires more power in the delta function than the other physical models.

Three models have an explicitly defined PAH fraction. The DL07 and HD21 models define the parameter q_PAH as the fraction of dust mass contained in carbonaceous grains with less than 10³ carbon atoms, roughly less than 1 nm in size. We use the f_PAH parameter from the MC11 model to estimate a PAH fraction, i.e., the mass of grains with sizes from 0.35 to 1.2 nm, as defined by the fiducial parameters. We refer to the mid-IR emission feature carriers in THEMIS as HACs. We use the definition in Lianou et al. (2019) to compute a fraction of HACs from THEMIS results, which can compare to the PAHs in other dust models: they found that this fraction of HACs corresponds to grains between 0.7 and 1.5 nm of the sCM20 component. It is important to remember that the strict definition of the PAH/HAC fraction is different in each model, but its purpose—fitting the mid-IR emission features—remains similar.

We investigate the variations of the surface density of the carriers, Σ_PAH and Σ_HAC, instead of their abundances. In the top panel of Figure 9, we show the radial profiles of Σ_{{PAH; HAC}}. Although the absolute values of the surface density of PAHs (HACs) differ by a factor of up to ∼3.5 (similar to dust masses), their gradients are very similar. This behavior shows that the grain populations that are held responsible for the mid-IR features in each model follow comparable distributions. In these models, their contributions to the total dust mass vary significantly but all prove to be a good fit to the mid-IR bands (see also Figure 5). This is also exemplified by the normalized ratio maps in Figure 9 (bottom panel). The dark colors in the outermost pixels of the HD21 model are due to the best q_PAH fit consistent with 0%. To visualize the variations between models, we normalize each parameter map to their mean value (as shown in the color-bar labels). We are thus able to compare the spatial variations of the maps and avoid the offsets due to the definition differences of PAHs or HACs.

Zoom In Zoom Out Reset image size

The q_PAH map in Aniano et al. (2020, using the DL07 model) shows similar features to ours. A large portion of the disk of M101 has a rather constant distribution of q_PAH, with conspicuous drops in H ii regions. In their study using the Desert et al. (1990) dust model, Relaño et al. (2020) found a flat radial profile of the small-to-large grain mass ratio, up to 0.8 r₂₅ (∼91). Our maps of the fraction of PAHs, or HACs, present a somewhat flat distribution (variations less than 1%) out to ∼0.3 r₂₅ (∼34; see Appendix B) and a steep change farther out.

To investigate in more detail the ability of each physical model to reproduce the PAH features, we perform a fit on an integrated SED and compare the results to measurements from the Infrared Spectrograph (IRS; on board Spitzer) in that same region (J. D. T. Smith 2021, private communication). Figure 10 shows a zoom on the mid-IR part of the models and the results of the fits to the integrated SED. From that fit, it appears that all models are able to generally reproduce the mid-IR features, with their different parameters, but we can notice a few differences between models.

Zoom In Zoom Out Reset image size

From the residuals (bottom panel), we see that the models perform similarly at 5.8, 8.0, and 12 μm. At 22 and 24 μm, the offset between the measurements means the models tend to split the difference and sit between the points. We note that all models appear to overestimate the continuum around 7 μm, in between the 6.2 and 7.7 μm PAH features.

There are nonetheless a couple of noticeable differences between each model. For instance, the HD21 model shows a higher continuum at 10 μm than the other three models, despite a similar continuum at 20 μm. This rules out the higher ${U}_{\min }$ found in the HD21 model as the reason for the higher flux at 10 μm. Rather, in this model, the emission at 10 μm is strongly dependent on the amount of nanosilicates used to account for the lack of correlation between PAH emission and anomalous microwave emission (Hensley & Draine 2017). The ratio between the flux at 20 and 10 μm is ∼1.7 for the HD21 models and 2 or above for the other three models. On the other hand, THEMIS shows no emission feature around 17 μm, while the other models do (although it has no impact on this particular fit). Around 7 μm, all models show a higher flux than the one seen in the IRS spectrum.

It is notable that, using only photometric bands from WISE and Spitzer /IRAC in the fit, all models reasonably reproduce well the mid-IR emission features, despite having different values of the PAH (or HAC) fraction and different definitions of the carriers. However, the comparison to spectroscopic measurements shows that there are still differences between models.

The monochromatic dust emission in far-IR wavelengths has often been used as a mass tracer of the ISM (e.g., Eales et al. 2012; Groves et al. 2015; Berta et al. 2016; Scoville et al. 2017; Aniano et al. 2020; Baes et al. 2020). In Figure 11, we plot the emission of M101 at 500 μm as a function of the fitted Σ_dust for each model (pixels above the 3σ detection threshold), color-coded by the minimum radiation field ${U}_{\min }$ or dust temperature T_dust.

Zoom In Zoom Out Reset image size

In all cases, we can see two distinct relations as the SPIRE 500 emission increases. The majority of the fitted pixels show a linear scaling between the emission at 500 μm and Σ_dust, while in some specific regions of the galaxy, all models prefer a higher radiation field (or temperature) and a lower dust surface density. We provide the scaling relations between the emission at 500 μm in MJy sr⁻¹, and the dust surface density in M_⊙ pc⁻², for each model. We measure the 5th and 95th percentiles of the data points above the 3σ detection threshold (to keep the bulk of the distribution only). We fit a linear slope to these points: ²⁴

$$\begin{eqnarray}\,\begin{array}{rcl}{{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{DL}07} & = & 1.21\pm 0.01\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.38\pm 0.005\\ {{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{MC}11} & = & 1.25\pm 0.02\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.39\pm 0.008\\ {{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{THEMIS}} & = & 1.32\pm 0.02\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.56\pm 0.008\\ {{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{HD}21} & = & 1.21\pm 0.01\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.50\pm 0.006\\ {{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{SE}} & = & 1.61\pm 0.03\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.80\pm 0.02\\ {{\mathrm{log}}_{10}({{\rm{\Sigma }}}_{d})}^{\mathrm{BE}} & = & 1.08\pm 0.08\times {\mathrm{log}}_{10}({I}_{\nu }^{500\,\mu {\rm{m}}})\\ & & -1.78\pm 0.005.\end{array}\,\end{eqnarray} \tag{ 15 }$$

To identify the pixels that "branch out" from the bulk, we select any pixel that falls below one standard deviation from the fit (dashed lines). In each panel, we show the spatial location of these pixels. It becomes clear that the regions that need a higher ${U}_{\min }$ (or T_dust) are H ii regions and in the outskirts of the galaxy. These pixels can account for between 4% and 11% of the pixels above the 3σ detection. The branching-out from the main relation is likely the consequence of the fact that the dust in these H ii regions is significantly hotter than average (as shown by the enhanced $\overline{U}$ in all models in Figure 8). Lower dust-to-gas ratios in the galaxy outskirts may also contribute to this trend.

In Aniano et al. (2020), the authors found that this relation is well represented by a power-law scaling in the KINGFISH sample, with a slope of 0.942, which is lower than our measured value of 1.2 for M101 with the DL07 model. A linear slope would be expected if the dust temperature and optical properties were uniform throughout the galaxy, leading to a constant ${I}_{\nu }^{500}/{{\rm{\Sigma }}}_{\mathrm{dust}}$ ratio. The spatial distribution of temperatures throughout each individual galaxy leads to distinct slopes, and for M101, we find that in regions of higher dust surface density, the dust is also warmer, leading to more 500 μm emission (this can be seen in the change in the color table in Figure 11 at the highest Σ_dust). The branch of H ii region points we see represents only a small fraction of the total data and may not be evident on the Aniano et al. (2020) plot.

The calibration of dust models involves a constraint on elements locked in grains (see Section 2.3). This step relies on depletion measurements, which characterize the distribution of heavy elements between the gas and solid phases. The final amount of elements allowed in dust grains varies between different physical dust models. The final dust masses derived by each model vary as well, as discussed in Sections 5.2 and 5.3.

A way to assess the performance of dust models is to verify that the required dust mass does not exceed the available heavy element mass, as constrained by metallicity measurements (e.g., Gordon et al. 2014; Chiang et al. 2018). We perform this test in M101 because its metallicity gradient has been thoroughly characterized (e.g., Zaritsky et al. 1994; Moustakas et al. 2010; Croxall et al. 2016; Berg et al. 2020; Skillman et al. 2020).

We estimate the dust mass surface density upper limit by assuming all available metals are in dust and calculating the metal mass surface density from the metallicity gradient and observed gas mass surface density:

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{\mathrm{dust}}^{\max }}{{M}_{\mathrm{gas}}}=\displaystyle \frac{{M}_{{\rm{Z}}}}{{M}_{\mathrm{gas}}},\end{eqnarray} \tag{ 16 }$$

where M_gas and M_Z are determined as follows.

The gas surface density is the sum of H i and H₂ surface densities including a correction for the mass of He (we ignore the ionized gas contribution). The latter is built from CO emission, assuming two prescriptions for the CO-to-H₂ conversion (see Section 2). We include the MW α_CO prescription as it is widely used (Equation (1)). We also choose the α_CO prescription from Bolatto et al. (2013), which takes into account environmental variations of α_CO with metallicity and surface density (Equation (2)). We emphasize that the Σ_dust upper limit in this section is dependent on the choice of α_CO to derive a gas surface density. This is particularly true in the central region of M101, where H₂ dominates (e.g., Schruba et al. 2011; Vílchez et al. 2019) and where the two α_CO differ the most. Note also that this result differs from that of Chiang et al. (2018). This is expected as the α_CO conversion factor in their study (from Sandstrom et al. 2013) is lower than the ones used in this study, which leads to a lower upper limit.

We use the $12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})$ radial profile from Berg et al. (2020) and convert it to metallicity through

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{{\rm{Z}}}}{{M}_{\mathrm{gas}}}=\displaystyle \frac{\tfrac{{m}_{{\rm{O}}}}{{m}_{{\rm{H}}}}\ {10}^{\left(12+{\mathrm{log}}_{10}({\rm{O}}/{\rm{H}})\right)-12}}{1.36\tfrac{{M}_{{\rm{O}}}}{{M}_{Z}}},\end{eqnarray} \tag{ 17 }$$

with m_O and m_H the atomic masses of oxygen and hydrogen, respectively, and the oxygen-to-metals mass ratio M_O/M_Z = 0.445 (Asplund et al. 2009).

The top panels in Figure 12 show the radial profile of Σ_dust for each model and the ${{\rm{\Sigma }}}_{\mathrm{dust}}^{\max }$ upper limits yielded by the two α_CO prescriptions used here (black lines, dotted and dashed–dotted). We see on the top-left panel that the DL07 and MC11 models are above both ${{\rm{\Sigma }}}_{\mathrm{dust}}^{\max }$ upper limits in almost all of the significant pixels. THEMIS and the HD21 models are fairly in line with the Bolatto et al. (2013) upper limit, with the conservative assumption that 100% metals are in dust grains. These behaviors suggest that the dust emissivity in all physical models is too low, which leads to requiring too much dust. We note that the large dust masses are likely due to the opacity calibration rather than a wrong fit in the far-IR bands: in Figure 5, all physical models show a reasonable fit at 160 μm, much closer to the IR peak than the 500 μm band. We are confident that the IR peak is correctly recovered and that the high dust masses are not due to the submillimeter excess. Based on the reasonable quality of the fits, we believe that the excessive mass is likely due to the opacity calibration.

Zoom In Zoom Out Reset image size

Despite sharing a common calibration approach, the details of the opacity calibration in the dust models used in this study vary in small but significant ways. While all models were calibrated to MW diffuse emission, they did not use exactly the same high-latitude cirrus spectrum. In addition, the M_dust/M_H adopted for the MW diffuse ISM by the physical dust models varies. Additionally, the radiation field that best reproduces the MW diffuse emission, U^MW, differs slightly from one model to the next. Because of the relationship between dust temperature and radiation field $(U\propto {T}_{\mathrm{dust}}^{4+\beta })$ and dust temperature and luminosity, even a slight difference in the assumed radiation field may lead to a significant change in the model's calibrated dust opacity.

To investigate calibration discrepancies, we renormalize each of the dust models via a fit to a common MW diffuse emission spectrum using the same abundance constraints. We use the MW SED described in Gordon et al. (2014), which we previously used to calibrate the κ_ν of the modified blackbody models (Chiang et al. 2018). This SED is the same as that used in Compiègne et al. (2011), a combination of DIRBE and FIRAS measurements (e.g., Boulanger et al. 1996; Arendt et al. 1998). ²⁵ We do not use the ionized gas correction because depletion measurements do not correct for it, and instead use a correction factor of 0.97 for molecular gas only (Compiègne et al. 2011). We integrate the SED in the PACS 100, PACS 160, SPIRE 250, SPIRE 350, and SPIRE 500 bands, so all models use the same wavelength coverage. We use 2.5% uncorrelated and 5% correlated errors to account for FIRAS and DIRBE uncertainties (Gordon et al. 2014). The M_dust/M_H ratio set in the normalization is 1/150, as suggested by depletion studies (F_* = 0.36 from Jenkins 2009, see also Gordon et al. 2014).

To perform the renormalization using the MW SED, we use the same fitting technique as previously described with the following choices: (1) we do not use the combination of radiation fields nor the stellar component (i.e., Ω_* = γ = 0). (2) We allow the minimum radiation field ${U}_{\min }$ and the total dust surface density Σ_dust to vary in each physical model. (3) We keep the relative ratios between grain populations fixed and do not vary them independently (e.g., for each model, we use the total spectra solid lines "DL07," "HD21," "MC11," and "THEMIS" in Figure 2).

The fits yield renormalization factors that correct all physical models from their respective assumptions to a dust-to-H ratio of 1/150. These corrections range from 1.5 for THEMIS, to 3 for the DL07 model (see Appendix A). With this normalization, we are able to meet the metallicity constraints. The top-right panel of Figure 12 shows the radial profiles with the correction factors applied to the surface densities of the physical dust models. The renormalization brings the models to lower dust surface densities that agree with the upper limit based on the metal content. It is interesting to note that the renormalized models now show three distinct behaviors in the dust mass surface density radial profile: DL07, HD21, and the broken power-law emissivity modified blackbody yield very similar results; THEMIS and MC11 are very similar to each other and offset by a factor of ∼2 from the first group, and the simple power-law emissivity modified blackbody has a steeper increase that still puts it above the abundance constraints even though it is similarly normalized to the MW cirrus spectrum. The two different behaviors, for DL07 and HD21, and for THEMIS and MC11 are likely due to the difference in the best ${U}_{\min }^{\mathrm{MW}}$ to fit the MW SED, and their initial spectrum used for calibration.

In the bottom panels of Figure 12, we show the radial profile of the dust-to-gas ratios (DGR) for each model, using the ${\alpha }_{\mathrm{CO}}^{\mathrm{BWL}13}$ conversion factor. The bottom left shows the DGR with the derived dust surface densities, and the bottom-right panels are after renormalization. The thick gray line shows the upper limit of the DGR using the metallicity gradient from Berg et al. (2020) and Equation (17), with ${\alpha }_{\mathrm{CO}}^{\mathrm{BWL}13}$. We find the same abrupt change in the DGR as Vílchez et al. (2019) around 0.5 r₂₅ , although with slightly lower values, consistent with the higher dust surface densities in this study.

DL07 was calibrated using the following constraints. The extinction is described in Weingartner & Draine (2001b) and uses the Fitzpatrick (1999) extinction curves with a normalization of N_H/E(B − V) = 5.8 × 10²¹ H cm⁻² or A_V /N_H = 5.3 × 10⁻²² cm². The high-latitude cirrus emission per H observed by DIRBE (Diffuse Infrared Background Experiment) and FIRAS (Far Infrared Absolute Spectrophotometer; Arendt et al. 1998; Finkbeiner et al. 1999) is used as a reference for the far-IR emission, complemented by mid- and near-IR emission from IRTS (Infrared Telescope in Space; Onaka et al. 1996; Tanaka et al. 1996). Weingartner & Draine (2001b) adopt solar abundances from Grevesse & Sauval (1998), assuming 30% of carbon is in the gas phase. They assume all silicon is depleted and has abundance equal to the solar value. DL07 uses M_dust/M_H = 1.0 × 10⁻². The radiation field used in the Draine & Li (2007) model is based on Mathis et al. (1983).

MC11 was calibrated using the following constraints. Extinction constraints were taken from Mathis (1990) and Fitzpatrick (1999), including the R_V = 3.1 extinction curve in the UV-visible and a normalization of N_H/E(B − V) = 5.8 × 10²¹ H cm⁻². At λ > 25 μm, MC11 use the MW cirrus emission per H observed by COBE-DIRBE and WMAP (integrated in the Herschel and Planck/HFI bands; see MC11). At λ ≤ 25 μm, a compilation of mid-IR observations of high-latitude MW cirrus is used (combining measurements from AROME, DIRBE, and ISOCAM; we refer to reader to the Compiègne et al. paper for details). They scale the emission SED by 0.77 to account for ionized and molecular gas not accounted for in the H column. The allowed dust-phase abundances for C, O, and other dust components come from the difference between solar (or F/G star) abundances and the observed gas-phase abundances. In total, M_dust/M_H = 1.02 × 10⁻². MC11 assumes the Mathis et al. (1983, D_G = 10 kpc) solar neighborhood radiation field to heat the dust grains.

THEMIS was calibrated using the same constraints as Compiègne et al. (2011) presented in the previous section, with the addition of the far-IR-to-extinction relation τ₂₅₀/E(B − V) = 5.8 × 10⁻⁴ (Planck Collaboration et al. 2011). In THEMIS, M_dust/M_H = 7.4 × 10⁻³.

The full set of observational constraints used to develop the model is described in Hensley & Draine (2021). In brief, the extinction curve is primarily a synthesis of those of Fitzpatrick et al. (2019) in the UV and optical, Schlafly et al. (2016) in the optical and near-infrared, and Hensley & Draine (2020) in the mid-infrared. The normalization ${N}_{{\rm{H}}}/E\left(B-V\right)=8.8\,\times {10}^{21}\,{\mathrm{cm}}^{-2}\,{\mathrm{mag}}^{-1}$ is used to normalize extinction to the hydrogen column (Lenz et al. 2017). The infrared emission in both total intensity and polarization are based on the analyses presented in Planck Collaboration Int. XVII (2014, 2015), and Planck Collaboration XI (2020), including the normalization to N_H. The solid phase interstellar abundances are re-determined in Hensley & Draine (2021) using a set of solar abundances (Asplund et al. 2009; Scott et al. 2015a, 2015b), a measurement of Galactic chemical enrichment from solar twin studies (Bedell et al. 2018), and determination of the gas-phase abundances from absorption spectroscopy (Jenkins 2009). M_dust/M_H = 1.0 × 10⁻² in this model. HD21 assumes the same radiation field as the DL07 model to heat the dust grains, with updates from Draine (2011).