The Resolved Distributions of Dust Mass and Temperature in Local Group Galaxies

TL;DR

This study maps dust temperature and mass distribution in four Local Group galaxies using Herschel data, revealing how resolution affects derived dust properties and mass estimates, especially in star-forming regions.

Contribution

It provides detailed dust temperature and mass surface density maps at high resolution and analyzes how resolution impacts dust property measurements across different galaxy types.

Findings

LMC has the highest mass-weighted average dust temperature.

Coarser resolution leads to systematically higher temperature estimates.

Total dust mass is underestimated at lower resolutions, especially in clumpy galaxies.

Abstract

We utilize archival far-infrared maps from the Herschel Space Observatory in four Local Group galaxies (Small and Large Magellanic Clouds, M31, and M33). We model their Spectral Energy Distribution (SED) from 100 to 500 $\mu$m using a single-temperature modified blackbody emission with a fixed emissivity index of $\beta = 1.8$. From the best-fit model, we derive the dust temperature, $T_{\rm d}$, and the dust mass surface density, $\Sigma_{\rm d}$, at 13 parsec resolution for SMC and LMC, and at 167 parsec resolution for all targets. This measurement allows us to build the distribution of dust mass and luminosity as functions of dust temperature and mass surface density. We compare those distribution functions among galaxies and between regions in a galaxy. We find that LMC has the highest mass-weighted average $T_{\rm d}$, while M31 and M33 have the lowest mass-weighted average $T_{\rm d}$. Within a galaxy, star forming regions have higher $T_{\rm d}$ and $\Sigma_{\rm d}$ relative to the overall distribution function, due to more intense heating by young stars and higher gas mass surface density. When we degrade the resolutions to mimic distant galaxies, the mass-weighted mean temperature gets warmer as the resolution gets coarser, meaning the temperature derived from unresolved observation is systematically higher than that in highly resolved observation. As an implication, the total dust mass is lower (underestimated) in coarser resolutions. This resolution-dependent effect is more prominent in clumpy star-forming galaxies (SMC, LMC, and M33), and less prominent in more quiescent massive spiral (M31).

Tables

Table 1. Table 1: Descriptions of symbols used in this paper

Symbols	Descriptions	Measured for	References
Temperature [K]
$T_d$	Expectation value of the dust equilibrium temperature at a given map resolution	a pixel	Equation 8
$T_{d,L}$	Luminosity-weighted of $T_d$ at a given map resolution	a pixel	Equation 14
${⟨T_d⟩}_M$	Mass-weighted mean of $T_d$ from a given map resolution	whole galaxy	Equation 16
${⟨T_d⟩}_L$	Luminosity-weighted mean of $T_d$ from a given map resolution	whole galaxy	Equation 17
$T_{\mathrm{peak}}$	Dust temperature that corresponds to the peak of mass or luminosity distribution	distribution function	Section 4.1
$T_{\mathrm{med}}$	The median value of dust temperature in the mass or luminosity distribution	distribution function	Section 4.1
${𝒲}_T$	The 16th-to-84th percentile range of temperature in mass or luminosity distribution	distribution function	Section 4.1
Mass [$M_{\odot }$] and Mass Surface Density [$M_{\odot }$pc-2]
${\Sigma }_d$	Expectation value of the dust mass surface density at a given map resolution	a pixel	Equation 9
${\Sigma }_{d,L}$	Luminosity-weighted of ${\Sigma }_d$ at a given map resolution	a pixel	Equation 15
$M_{d,\mathrm{Tot}}$	Total mass from integrated SED using single- or multi-temperature model	whole galaxy	Section 6.2
${\Sigma }_{\mathrm{peak}}$	Dust mass surface density that corresponds to the peak of mass or luminosity distribution	distribution function	Section 4.2
${\Sigma }_{\mathrm{med}}$	The median value of dust mass surface density in the mass or luminosity distribution	distribution function	Section 4.2
${𝒲}_{\Sigma }$	The 16th-to-84th percentile range of mass surface density in mass or luminosity distribution	distribution function	Section 4.2
Interstellar Radiation Field
$U$	Interstellar radiation field strength (energy density) at a given resolution relative to Solar neighborhood value	a pixel	Equation 10
$U_{\min }$	The minimum value of $U$ used in multi-temperature modeling	model	Equation 18
$U_{\max }$	The maximum value of $U$ used in multi-temperature modeling	model	Equation 18
$U_{\mathrm{peak}}$	Interstellar radiation field strength that corresponds to the peak of mass or luminosity distribution	distribution function	Section 4.1
$U_{\mathrm{med}}$	The median value of $U$ in the mass or luminosity distribution	distribution function	Section 4.1
${𝒲}_U$	The 16th-to-84th percentile range of $U$ in mass or luminosity distribution	distribution function	Section 4.1
Others
${𝒥}_{\mathrm{IR}}$	A quantity proportional to the IR equilibrium luminosity at a given map resolution	a pixel	Equation 11
$\alpha$	The slope of the power-law distribution of dust mass heated by $U$	model	Dale et al. (2001)
$\beta$	The dust emissivity index (power law exponent in the dependency of $\kappa$ as a function of wavelength)	model	Equation 4
$\gamma$	The fraction of dust mass heated by a power law distribution field, $U^{-\alpha }$, between $U_{\min }$ and $U_{\max }$	model	Draine & Li (2007)
$\kappa$	Dust absorption cross-section	model	Equation 3
$ℒ$	Likelihood of a given model	model	Equation 6
$r_s$	Spearman rank correlation coefficient	data	Spearman (1904)

Table 2. Table 2: The median value and the median absolute deviation (MAD) of the 500 µm map.

	SMC	LMC	M31	M33
Median [MJy sr-1]	0.012	0.187	1.663	0.135
MAD [MJy sr-1]	0.215	0.493	0.212	0.210

Table 3. Table 3: Parameters of the distribution function at 13 pc and 167 pc resolutions with β = 1.8 𝛽 1.8 \beta=1.8

Galaxies	Resolutions
		Radiation Field			Dust Equilibrium Temperature			Dust Mass Surface Density
		log${}_{10}{}^{}U_{\mathrm{peak}}^{}$	log${}_{10}{}^{}U_{\mathrm{med}}^{}$	Width [dex]	$T_{\mathrm{peak}}$ [K]	$T_{\mathrm{med}}$ [K]	Width [K]	log${}_{10}{}^{}\Sigma _{\mathrm{peak}}^{}$	log${}_{10}{}^{}\Sigma _{\mathrm{med}}^{}$	Width [dex]
Distribution of Mass
SMC	13 pc	$0.03$	$0.01$	$0.66$	$18.18$	$18.08$	$4.74$	$-1.36$	$-1.33$	$0.67$
LMC		$0.27$	$0.23$	$0.59$	$20.00$	$19.73$	$4.67$	$-1.27$	$-1.24$	$0.78$
SMC	167 pc	$0.16$	$0.15$	$0.43$	$19.14$	$19.10$	$3.28$	$-1.74$	$-1.49$	$0.66$
LMC		$0.30$	$0.28$	$0.48$	$20.24$	$20.15$	$3.92$	$-1.33$	$-1.34$	$0.68$
M31		$-0.01$	$-0.14$	$0.60$	$17.89$	$17.00$	$3.97$	$-1.69$	$-1.65$	$0.58$
M33		$-0.52$	$-0.33$	$1.04$	$14.61$	$15.76$	$6.72$	$-1.31$	$-1.36$	$0.48$
Distribution of Luminosity
SMC	13 pc	$0.18$	$0.26$	$0.68$	$19.30$	$19.94$	$5.45$	$-1.15$	$-1.26$	$0.90$
LMC		$0.34$	$0.47$	$0.74$	$20.56$	$21.66$	$6.52$	$-1.11$	$-1.07$	$0.95$
SMC	167 pc	$0.18$	$0.23$	$0.37$	$19.30$	$19.76$	$2.93$	$-1.35$	$-1.38$	$0.69$
LMC		$0.35$	$0.44$	$0.60$	$20.64$	$21.42$	$5.23$	$-1.28$	$-1.21$	$0.74$
M31		$0.04$	$0.03$	$0.57$	$18.25$	$18.21$	$4.12$	$-1.43$	$-1.63$	$0.67$
M33		$0.44$	$0.27$	$0.82$	$21.39$	$20.00$	$6.23$	$-1.19$	$-1.24$	$0.45$

Table 4. Table 4: Summary of measurements for the whole galaxy.

Galaxies	Flux Density [${10}^3$ Jy]					Dust Equilibrium Temperature [K]			Total Dust Mass [${10}^5{M}_{\odot }$]
	100 µm	160 µm	250 µm	350 µm	500 µm	Single $T_d$	Multi $T_d$	Resolved	Single $T_d$	Multi $T_d$	Resolved
SMC	$14.07\pm 1.69$	$18.73\pm 2.25$	$9.64\pm 0.92$	$5.61\pm 0.53$	$2.65\pm 0.25$	${19.6}_{1.2}^{0.7}$	${17.6}_{-1.0}^{-1.2}$	$18.2\pm 1.7$	${1.6}_{0.3}^{0.3}$	${3.9}_{1.1}^{1.3}$	$1.9\pm 0.7$
LMC	$222.68\pm 26.72$	$264.83\pm 31.78$	$138.42\pm 13.15$	$71.19\pm 6.76$	$30.13\pm 2.86$	${21.3}_{1.1}^{0.6}$	${21.2}_{-1.8}^{-1.7}$	$20.0\pm 1.3$	${11.7}_{1.9}^{1.5}$	${15.9}_{5.2}^{5.9}$	$13.1\pm 3.5$
M31	$3.59\pm 0.43$	$8.71\pm 1.05$	$5.66\pm 0.54$	$3.05\pm 0.29$	$1.32\pm 0.13$	${17.4}_{0.7}^{0.1}$	${18.1}_{-0.9}^{-0.7}$	$16.9\pm 1.1$	${195.7}_{26.2}^{17.9}$	${192.1}_{22.0}^{8.7}$	$195.0\pm 46.5$
M33	$1.46\pm 0.18$	$2.14\pm 0.26$	$1.27\pm 0.12$	$0.71\pm 0.07$	$0.33\pm 0.03$	${19.2}_{1.0}^{0.5}$	${18.4}_{-1.6}^{-1.4}$	$16.3\pm 1.3$	${41.0}_{7.1}^{5.9}$	${73.2}_{21.2}^{25.5}$	$61.7\pm 23.3$

Table 5. Table 5: The expected values of integrated SED multi-temperature fits.

Galaxies	$\alpha$	log${}_{10}{}^{}U_{\min }^{}$	$\gamma$
		[dex]	[%]
SMC	${1.61}_{-0.45}^{+0.19}$	$-{0.51}_{-0.32}^{+0.12}$	${20.73}_{-16.16}^{+26.80}$
LMC	${2.15}_{-0.95}^{+1.01}$	${0.12}_{-0.38}^{+0.18}$	${8.78}_{-7.50}^{+29.36}$
M31	${3.31}_{-1.29}^{+1.08}$	$-{0.11}_{-0.17}^{+0.03}$	${1.51}_{-1.36}^{+8.88}$
M33	${1.77}_{-0.60}^{+0.20}$	$-{0.32}_{-0.33}^{+0.13}$	${14.27}_{-11.49}^{+29.65}$

Equations

$$I_{\nu }(\lambda )=B_{\nu }(T_{\mathrm{d}},\lambda )[1-e^{-\tau }],$$

$$I_{\nu }(\lambda )\approx \tau (\lambda )B_{\nu }(T_{\mathrm{d}},\lambda ).$$

$$\tau (\lambda )=\kappa (\lambda ){\Sigma }_{\mathrm{d}},$$

$$\kappa (\lambda )=\kappa ({\lambda }_0){\left(\frac{\lambda }{{\lambda }_0}\right)}^{-\beta },$$

$$I_{\nu }(\lambda )=\kappa ({\lambda }_0){\Sigma }_{\mathrm{d}}{\left(\frac{\lambda }{{\lambda }_0}\right)}^{-\beta }{B}_{\nu }(T_{\mathrm{d}},\lambda ).$$

$$\mathcal{L}=\exp (-{\chi }^2/2),$$

$${\chi }^2=\Delta I^T{\mathcal{C}}^{-1}\Delta I,$$

$$T_{\mathrm{d}}=\frac{{\sum }_i{\mathcal{L}}_i{T}_{\mathrm{d},i}}{{\sum }_i{\mathcal{L}}_i},$$

$${\log }_{10}{\Sigma }_{\mathrm{d}}=\frac{{\sum }_i{\mathcal{L}}_i{\log }_{10}{\Sigma }_{\mathrm{d},i}}{{\sum }_i{\mathcal{L}}_i}.$$

$${\log }_{10}U=(4+\beta ){\log }_{10}\left(\frac{T_{\mathrm{d}}}{18\mathrm{K}}\right).$$

$${\mathcal{J}}_{\mathrm{IR}}\equiv {\Sigma }_{\mathrm{d}}{T}_{\mathrm{d}}^{4+\beta }.$$

$$\frac{d{m}_{\mathrm{d}}}{d{\log }_{10}U}(U)=\frac{1}{M_{\mathrm{Tot}}}\frac{d{M}_d}{d{\log }_{10}U}(U).$$

$$\frac{d{m}_{\mathrm{d}}}{d{\log }_{10}{T}_{\mathrm{d}}}=(4+\beta )\frac{d{m}_{\mathrm{d}}}{d{\log }_{10}U}.$$

$$T_{d,L}=\frac{{\mathcal{J}}_{\mathrm{IR}}{T}_d\ast \mathcal{G}}{{\mathcal{J}}_{\mathrm{IR}}\ast \mathcal{G}},$$

$${\Sigma }_{d,L}=\frac{{\mathcal{J}}_{\mathrm{IR}}{\Sigma }_d\ast \mathcal{G}}{{\mathcal{J}}_{\mathrm{IR}}\ast \mathcal{G}}.$$

$$⟨T_{\mathrm{d}}{⟩}_M=\frac{{\sum }_i{\Sigma }_{d,L,i}{T}_{\mathrm{d},L,i}}{{\sum }_i{\Sigma }_{d,L,i}},$$

$$⟨T_{\mathrm{d}}{⟩}_L=\frac{{\sum }_i{\mathcal{J}}_{\mathrm{IR},\mathrm{i}}{T}_{\mathrm{d},L,i}}{{\sum }_i{\mathcal{J}}_{\mathrm{IR},\mathrm{i}}}.$$

$$\frac{d{m}_d}{dU}=(1-\gamma )\delta (U-U_{\min })+\gamma \frac{(\alpha -1)}{U_{\min }^{1-\alpha }-U_{\max }^{1-\alpha }}{U}^{-\alpha },$$

$$⟨U{⟩}_M=(1-\gamma )U_{\min }+\gamma \frac{U_{\min }\ln (U_{\max }/U_{\min })}{1-U_{\min }/U_{\max }}.$$

Full text

ISM: dust — infrared: ISM — galaxies: ISM — galaxies: Local Group

††facilities: Herschel(PACS and SPIRE), Spitzer(MIPS).††software: CASA (v4.7.2; McMullin et al., 2007).

1 Introduction

Most of the dust mass in galaxies resides in grains that are in thermal equilibrium with the interstellar radiation field (ISRF; Draine, 2003). The strength of the interstellar radiation field is often denoted as $U$, which is the energy density of starlight relative to that measured by Mathis et al. (1983) for the Solar neighborhood. For a single population of dust in equilibrium with a radiation field, $U$, the dust temperature, $T_{\rm d}$, depends on the radiation field in a simple way via $U\propto T_{\rm d}^{4+\beta}$, where $\beta$ is the dust emissivity index. That temperature-ISRF relation assumes the dust cross-section per unit mass, $\kappa$, depends on the wavelength as $\kappa\propto\lambda^{-\beta}$ (Draine, 2011).

Extragalactic observations necessarily convolve together many different environments and radiation fields due to the limited angular resolution of infrared telescopes. To account for this, starting from the work of Dale et al. (2001), many extragalactic studies have modelled the infrared spectral energy distribution (SED) using a combination of dust populations, each in equilibrium with a distinct radiation field (e.g., Dale & Helou, 2002; Draine et al., 2007; Galliano et al., 2011; Chastenet et al., 2017). In this model, a distribution function describes the amount of dust mass, $M_{d}$, heated by radiation fields, as $dM_{\rm d}/dU$, and is usually modeled as a power law, $U^{-\alpha}$, where $\alpha\sim 2$ (Draine et al., 2007).

Some part of this distribution of dust heating, $dM_{\rm d}/dU$, reflects geometry and radiative transfer on very small scales. For example, in the sub-parsec zone of influence of a young stellar population or the photon-dominated region at the edge of a molecular cloud, dust is exposed to a wide range of radiation field intensities (e.g., see Dale et al., 2001). These effects can be best studied in analytical radiative transfer models, highly resolved observations of Milky Way regions, or via their imprints on the integrated SED.

On larger scales, some part of the dust heating distribution will also arise from variations in the physical conditions across a galaxy (e.g., Gordon et al., 2014). The magnitude of these variations can be measured by mapping the resolved temperature and mass of dust across a galaxy. From such maps, we can measure how dust mass is distributed as a function of the illuminating radiation field. Also, we can estimate the effect of blurring together all the regions of a galaxy into a single measurement, which is equivalent to observing a more distant galaxy.

Due to their proximity, Local Group galaxies offer the best opportunity to study the resolved distributions of dust mass and dust temperature across the entire galaxy. In this paper, we use data from the Herschel Space Observatory (Pilbratt et al., 2010) to derive highly resolved maps ($13-167$ pc resolution) of dust mass surface density and temperature for four Local Group galaxies, the Small and Large Magellanic Clouds (SMC and LMC), M31, and M33. We use these maps to measure how dust mass is distributed as a function of dust temperature, which traces the average illuminating radiation field $U$ in an observational element or a pixel. We compare these distributions among our four targets, and we investigate how degrading the resolution of the data (equivalent to integrating over a range of dust mass and temperatures) would affect the inferred values of the dust mass and temperature at that degraded resolution. We also compare it with the dust temperature and dust total mass derived from integrated SED of unresolved object, that mimics the high-redshift studies (e.g., Magdis et al., 2010, 2012; Scoville et al., 2014; Genzel et al., 2015).

A number of previous papers have utilized the Herschel maps of dust emission in each of our targets galaxies individually. These include the studies of the SMC and LMC by Meixner et al. (2013); Gordon et al. (2014); Roman-Duval et al. (2014); Chastenet et al. (2017), the M31 work by Fritz et al. (2012); Groves et al. (2012); Smith et al. (2012); Draine et al. (2014), and the M33-focused investigations of Braine et al. (2010); Boquien et al. (2011); Xilouris et al. (2012). The new contributions of this paper include homogenizing the methodology and analysis for all targets, focusing on implications of the small scale structures ($\sim 10-100$ pc) for observations of more distant galaxies, and employing the new fitting code developed by Chiang et al. (2018), based on Gordon et al. (2014).

Throughout this paper, we use the five longest wavelengths of Herschel bands (100 µm, 160 µm, 250 µm, 350 µm, and 500 µm). We expect the infrared emission at these wavelengths mainly captures emission from relatively large grains in thermal equilibrium with the local radiation field. These grains tend to represent the dominant mass component, and focusing on this wavelength range allows us to employ a modified blackbody model to fit the SED, following Chiang et al. (2018) and Gordon et al. (2014) methodologies.

This paper is organized as follows. We describe the archival data in $\S$2 and our modified blackbody modeling in $\S$3. We present maps of dust mass surface density and temperature for all four targets in $\S$4. We also show the distribution of dust mass as functions of illuminating radiation field for each galaxy and locations within a galaxy in $\S$4. We show the correlation between dust mass surface density vs. dust temperature, IR color vs. dust temperature, and 500 µm intensity vs. dust mass surface density in $\S$5. We explore how the inferred dust temperature changes as a function of resolution (down to treating the galaxy as one unresolved source) in $\S$6. Finally, we summarize our findings in $\S$7.

Throughout this work, we adopt a distance of 62.1 kpc to the SMC (Graczyk et al., 2014), 50.2 kpc to the LMC (Klein, 2014), 744 kpc to M31 (Vilardell et al., 2010), and 840 kpc to M33 (Freedman et al., 1991). We also adopt inclination of $2^{\circ}.6$ for SMC (Subramanian & Subramaniam, 2012), $34^{\circ}.7$ for LMC (van der Marel & Cioni, 2001), $77^{\circ}.7$ for M31 (Corbelli et al., 2010), and $56^{\circ}$ for M33 (Paturel et al., 2003). We summarize the symbols used in this paper in Table 1.

2 Data

2.1 Herschel Far-infrared Maps

We utilize the far-infrared maps at $\lambda=100-500$ µm as observed by the Herschel Space Observatory (Pilbratt et al., 2010) using the PACS (Poglitsch et al., 2010) and SPIRE (Griffin et al., 2010) instruments. These maps were the data products of several key projects: HERITAGE used Herschel to observe the SMC and LMC (Meixner et al., 2013; Gordon et al., 2014; Roman-Duval et al., 2014), Groves et al. (2012) and Draine et al. (2014) present Herschel observations of M31 (see also Smith et al., 2012), and HerM33s (Braine et al., 2010; Boquien et al., 2011; Xilouris et al., 2012) used Herschel to observe M33.

Before any analysis, we convolve each map to convert the point spread function (PSF) from the original PSF of the instrument to a symmetric Gaussian PSF. We use the kernel that is appropriate for each PACS and SPIRE band (Aniano et al., 2011) to do this conversion. This Gaussian PSF allows easy convolution to any lower spatial resolution using coarser Gaussian kernels. After PSF conversion, these PACS and SPIRE maps have full width at half max (FWHM) resolutions of 15″ for 100 µm and 160 µm, 30″ for 250 µm and 350 µm, and 41″ for 500 µm. We refer to these values as the “native resolutions.”

2.2 Region of Interest

We identify a region of interest for each galaxy based on the SPIRE maps. We use the SPIRE maps to identify the region of interest because these maps have the best sensitivity to low dust column densities and cool dust. This region represents our estimate of the full spatial extent of the galaxy as detected by Herschel. We avoid this region of interest when fitting foreground and background emission, and we only fit the resolved SED within this region of interest.

To identify the region of interest, we first calculate the median value away from the galaxy and median absolute deviation (MAD) in the 500 µm map. Then, we define a threshold for significant emission equal to the median value plus $3\times$ the noise level as inferred from the median absolute deviation (MAD, Table 2). We consider the contiguous region of pixels with intensity above this threshold that is also contiguous with the main body of the galaxy to be our region of interest. We further dilate the region of interest by several beam widths after applying the threshold. This dilation includes faint emission around the galaxy and remove any holes in the region of interest.

These regions of interest do an excellent job of covering the whole body of each galaxy (shown as white contours in Figure 1). Given this, the details of this masking process have little impact on our overall results. We verified this by changing our adopted threshold to use a signal to noise cut of $2$ or $4$ and also changing the number of iterations to grow the mask. Modest changes only slightly affect the low-end tail of the distribution of dust mass surface density. The total dust mass in LMC, M31, and M33 only varies within few percent by using smaller or bigger mask. The region of interest definition has more effect on SMC, which is surrounded by extended, faint tidal features. As a result, changing the region of interest in the SMC can cause the total mass in our analysis to vary by $10-20\%$. As Figure 1 shows, our adopted region of interest does a good job of including the main body of the galaxy and the bright regions of the SMC’s eastern wing.

2.3 Foreground and Background Subtraction

Each Herschel image blends emission from the galaxy with foreground and background emission. These foreground and background contributions reflect a mixture of Milky Way cirrus, the cosmic infrared background and resolved background galaxies, and the adopted observing and imaging strategies which can result in zero-point offsets and gradients. We correct each map for this contamination in several steps.

First, following Bot et al. (2004), we use a Milky Way 21-cm map as a template to correct the SMC and LMC maps for foreground emission from the Galactic dust. Because of the wide sky coverage of these two galaxies, using Galactic 21-cm map represents the best way to capture foreground variations across those galaxies. We use the foregrounds calculated by Chastenet et al. (2017) using the Hi maps from Stanimirovic et al. (1999) and Staveley-Smith et al. (2003). We refer to Chastenet et al. (2017) for more details.

After the first step, for all galaxies and all bands, we estimate the combined foreground and backgrounds emission (hereafter just “background”) and subtract it from the map. To do so, we calculate the median intensity of the image outside the region of interest. We subtract this value from the whole image. We further refined this background estimate by fitting a plane outside the region of interest. We iterate this fit, and dropping pixels that deviated by more than $\pm 2\sigma$ from the fit, where $\sigma$ is the robustly estimated RMS noise. Exclusion of these outlier pixels is important to exclude pixels that may be contaminated by instrumental artifacts or non-background emission. Then, we subtract the fit plane from the whole image. Finally, we construct a histogram of intensities outside the region of interest and make one final adjustment to the zero point of the image, forcing it to coincide with the mode of the histogram. This final step usually represent a small adjustment (significantly less than $1\sigma$).

After this step, the zero point of all bands and all galaxies outside the region of interest appears reasonably consistent with zero intensity.

2.4 Convolution and Reprojection

We convolve the background-subtraced maps from their native angular resolutions to a set of common physical resolutions at the adopted distance of each galaxy. For the SMC and the LMC, we create maps with FWHM from 13 to 500 pc resolution. For M31 and M33 we create maps with FWHM from 167 to 2,000 pc resolution. In the following analysis, we work mostly at the finest physical scale, i.e. 13 pc for SMC and LMC, and 167 pc for all targets. We choose this scale to match the angular resolutions of SPIRE at 500 µm at the distance to the SMC and M33. Thus, the finest common resolution to study all four galaxies, while still including all SPIRE bands, is 167 pc.

After convolution, we reproject the images from all bands onto a common coordinate grid. We choose the pixel spacing for the common grid so that there are 2.5 pixels across each (FWHM) beam, i.e. we oversample at roughly the Nyquist sampling rate. This means the pixel size in parsec is also larger in coarser resolution. To carry out the reprojection and smoothing, we use the CASA tasks imsmooth and imregrid (McMullin et al., 2007).

3 Modeling

3.1 Modified Blackbody Model

The intensity of dust emission, $I_{\nu}$, for dust in equilibrium with the radiation field depends on the dust optical depth, $\tau(\lambda)$, and the dust equilibrium temperature, $T_{\rm d}$, via

[TABLE]

where $B_{\nu}(T_{\rm d},\lambda)$ is the Planck function for temperature $T_{\rm d}$ at a wavelength $\lambda$.

For our target galaxies and resolutions, dust emission at $\lambda\geq 100~{}\micron$ is almost always optically thin. In this case,

[TABLE]

The dust optical depth can be written as

[TABLE]

where $\kappa(\lambda$) is the dust absorption cross-section per unit mass at a wavelength $\lambda$, and $\Sigma_{\rm d}$ is the dust mass surface density.

Fitting the SED yields $\tau(\lambda)$, but for a given $\tau(\lambda)$, the values of $\kappa(\lambda)$ and $\Sigma_{\rm d}$ are degenerate. Breaking this degeneracy requires an independent calibration of $\kappa(\lambda)$. Such a calibration, in turn, requires that $\tau(\lambda)$ be measured in a location where $\Sigma_{\rm d}$ is known from independent measurements. The most common site for such calibrations is the local diffuse ISM in the Milky Way, where the spectroscopic measurements of depletion into dust yield an independent estimate of the dust mass surface density (e.g., Jenkins, 2009), while the SED has been measured by all-sky IR mapping missions.

We follow this approach, by adopting the calibration scheme of Chiang et al. (2018), based on Gordon et al. (2014). We calculate $\tau(\lambda)$ by fitting the Solar Neighborhood cirrus emission using our modified blackbody model. Then, we derive $\kappa(\lambda)$ from $\tau(\lambda)$ by adopting the dust abundance found by Jenkins (2009). We further assume that $\kappa(\lambda)$ remains constant across all regions in Local Group galaxies. Though variations in $\kappa$ do represent a systematic uncertainty, this procedure removes another systematic uncertainty by using the same dust model to estimate $\kappa$ and fit the data.

For dust, the absorption cross section per unit mass decreases with increasing wavelength. We follow standard practice and assumed it to be a power law (Draine, 2011) as

[TABLE]

where $\lambda_{0}$ is a reference wavelength and $\beta$ represents the dust emissivity index. In this paper, we assume a fixed $\beta=1.8$ (Dunne & Eales, 2001; Draine et al., 2007; Clements et al., 2010; Planck Collaboration et al., 2011; Scoville et al., 2014). According to Chiang et al. (2018), the variation in $\Sigma_{d}$ from different choices of $\beta$ values would be small since $\kappa(\lambda_{0})$ is calibrated accordingly.

Combining those assumptions above, we get

[TABLE]

We adopt $\lambda_{0}=160~{}\micron$ as our reference wavelength. For $\beta=1.8$, our fit to the Milky Way cirrus yields $\kappa(\lambda_{0})=18.7\pm 0.6$ cm2 g*-1*.

3.2 Fitting the Dust SED

For each line-of-sight, we estimate $\Sigma_{\rm d}$ and $T_{\rm d}$ by simultaneously fitting the intensities at 100 µm, 160 µm, 250 µm, 350 µm, and 500 µm using the model in Equation 5. We follow the algorithms presented in Chiang et al. (2018) and Gordon et al. (2014) to calculate the relative likelihood of the model, $\mathcal{L}$, given the observed SED. Here,

[TABLE]

and

[TABLE]

where $\Delta I$ indicates a vector recording the difference between the observed SED and the model SED. Before calculating $\Delta I$, we integrate the models (Equation 5) with the appropriate response function of that Herschel band.

$\mathcal{C}^{-1}$ refers to the the inverse of the covariance matrix. This term takes into account the inter-band covariances in the uncertainty calculation. T denotes matrix transpose operation. In the Gordon et al. (2014) approach (upon which Chiang et al. (2018) is based), this covariance matrix, $\mathcal{C}$, is calculated by summing the covariance matrix from the data outside the region of interest and the covariance matrix of the calibration uncertainty from the Herschel instruments (Balog et al., 2014; Bendo et al., 2017). The covariance of the data outside the region of interest are taken to indicate the covariance in the observational noise.

The fitting process calculates the likelihood, $\mathcal{L}$, for each model across a grid of $T_{\rm d}$ and $\Sigma_{\rm d}$. The grid for $T_{\rm d}$ spans from 5 to 50 K with increments of 0.5 K. We space the grid for $\Sigma_{\rm d}$ logarithmically, covering from $-4$ to $1$ dex in log${}_{10}\Sigma_{\rm d}$ with a step size of 0.025 dex.

From the grid of likelihoods, we calculate the expectation values for $T_{\rm d}$ and $\Sigma_{\rm d}$ via

[TABLE]

and

[TABLE]

Here, the sum, $\sum_{i}$, covers all cells in the grid. We adopt these expectation values as our estimate of the best-fit $\Sigma_{\rm d}$ and $T_{\rm d}$ at each pixel inside the region of interest. We correct $\Sigma_{d}$ for galaxy inclination, $i$, by multiplying $\Sigma_{d}$ with cos($i$).

We repeat this fit and the covariance estimation, at several resolutions. For this paper the key resolutions are $13$ pc (the finest common resolution achievable for both the SMC and LMC) and $167$ pc (the finest common resolution achievable for all four galaxies).

3.3 Interstellar Radiation Field Strength and Equilibrium Dust Luminosity

For dust in thermal equilibrium with the local radiation field and with a power-law mass absorption coefficient $\kappa$, $T_{\rm d}$ relates to the interstellar radiation field (ISRF), $U$, via

[TABLE]

The normalization takes $T_{\rm d}=18$ K to corresponds to the Solar Neighborhood radiation field, $U_{\odot}=1$ (Draine et al., 2014). This normalization agrees well with our fit to calibrate $\kappa$. In that case, we found $T_{\rm d}=18.29\pm 0.11$ K for the local cirrus. Assuming that $U_{\odot}$ heats local cirrus, then this supports equating $T_{\rm d}\approx 18$ K with $U\approx 1$.

The infrared luminosity emerging from grains that are in equilibrium with the local radiation field is the integral of Equation 2 over infrared wavelengths. For simplicity, we define a quantity proportional to the infrared luminosity per unit area as

[TABLE]

Here, $\mathcal{J}_{\rm IR}$ represents the infrared luminosity surface density of grains with mass surface density $\Sigma_{\rm d}$ at temperature $T_{\rm d}$. We refer to $\mathcal{J}_{\rm IR}$ as the “equilibrium luminosity surface density.”

Even though the infrared luminosity captures the integral of infrared SED under our best-fit model, it underestimates the true total infrared luminosity because we do not fit or consider emission from hot grains and/or small grains that are out of equilibrium with the local radiation field. Instead, it captures the total light emerging from the large grain population in equilibrium with the local ISRF. Contributions from these hot grains and small grains that are out of equilibrium only make a small difference in the total dust mass. However, they produce an important fraction of the total infrared luminosity emitted by all kind of dust.

3.4 Caveats

We adopt several simplifying assumptions that the reader should bear in mind. We assume a single $T_{\rm d}$ for each pixel. In reality, $T_{\rm d}$ may vary within any given resolution element, especially in star-forming regions, which will contain many local heating sources. Hence, our $T_{\rm d}$ measured at 13 or 167 pc resolution is a luminosity weighted mean of the sub-resolution $T_{\rm d}$ distribution. Our paper addresses the large-scale averaging that depends mainly on galaxy structure and the location of star-forming regions. The properties of the small scale distribution of $T_{\rm d}$ (Draine & Li, 2007) are beyond the scope of our work, and will be addressed in upcoming paper (Chastenet et al. in preparation). We refer the reader to Bianchi (2013) who performed a comparison between the single-temperature model and the full dust distribution model.

Variations in $\beta$ also affect the best-fit $T_{\rm d}$, with the sense that a higher value of $\beta$ gives a lower value of the best-fit $T_{\rm d}$ (Kelly et al., 2012). There are some measurements suggesting a variation of $\beta$ within and among Local Group galaxies (e.g., Smith et al., 2012; Gordon et al., 2014). We fix $\beta$ in the interests of simplicity and to avoid the ambiguities in the interpretation, as highlighted by Kelly et al. (2012), but this remains an important open topic.

For the dust mass, our measurements depend on the calibrated value of $\kappa(\lambda_{0})=18.7\pm 0.6$ cm2 g*-1* for $\lambda_{0}=160~{}\micron$ and $\beta=1.8$. This value relies on the assumption that the optical properties of dust grains throughout the Local Group galaxies is the same as in the Milky Way cirrus. If $\kappa(\lambda_{0})$ is underestimated, we will overestimate the dust mass, and vice-versa. In order to partially control for this uncertainty, we often normalize our results by the total dust mass in each galaxy. In the case that $\kappa(\lambda_{0})$ varies from galaxy-to-galaxy, but not between regions within a galaxy, our results will still be robust.

4 Maps and Distributions of Dust Mass, Luminosity, and Temperature

Figure 1 shows our best estimates of dust equilibrium temperature, $T_{\rm d}$, equilibrium IR luminosity surface density, $\mathcal{J_{\rm IR}}$, and dust mass surface density, $\Sigma_{\rm d}$, for each target.111These maps are available in FITS file at https://www.asc.ohio-state.edu/astronomy/dustmaps/ In the SMC, LMC, and M33, we find high dust temperatures in regions of active star formation, consistent with previous works (e.g., Xilouris et al., 2012; Gordon et al., 2014). In M31, the highest dust temperatures coincide with the low dust mass surface density in the inner, bulge-dominated part of the galaxy. Previous work has demonstrated that this hot dust in the inner part of M31 results from heating by the old stellar population (Groves et al., 2012; Draine et al., 2014; Viaene et al., 2014).

The dust mass surface density maps show the same features as gas mass surface density maps in these galaxies. In the SMC, the ISM material is concentrated along the bar, with an extension to the east into the wing, hosting the bright star-forming complex N83/N84. The LMC shows a prominent ridge of material south of 30 Doradus along the eastern edge of the galaxy and shells through the rest of the galaxy. The dust in M31 is concentrated into a series of rings that may be tightly wound spirals (e.g., Nieten et al., 2006; Gordon et al., 2006). And in M33, the dust mass maps show flocculent spiral structures.

Atomic hydrogen, Hi, makes up most of the neutral ISM in each of our targets. As a result, our $\Sigma_{\rm d}$ maps resemble Hi 21-cm maps of these galaxies (see Kim et al., 1998; Stanimirovic et al., 1999; Staveley-Smith et al., 2003; Braun et al., 2009; Braun, 2012; Koch et al., 2018). In detail, our maps should also reflect the presence of molecular gas and variations of the dust-to-gas ratio. Even accounting for dust-to-gas ratio variations, Figure 1 gives among the most uniform views of the ISM in Local Group galaxies up to date.

The luminosity maps ($\mathcal{J_{\rm IR}}\equiv\Sigma_{d}T_{d}^{5.8}$) show the brightest regions where most infrared light comes from. These bright regions are due to high $\Sigma_{d}$, high $T_{d}$, or both (excluding emission from hot, small grains). The star forming regions are still bright because they have both warm $T_{\rm d}$ and dense $\Sigma_{d}$. M31 shows distinctive feature, where the bulge (that shows a ‘hole’ in the $\Sigma_{d}$ map) is luminous because of its high $T_{\rm d}$.

Based on these maps, we calculate the distributions of dust mass and luminosity as functions of local illuminating radiation field, $U\propto T_{d}^{5.8}$ (in $\S$4.1), and local dust mass surface density, $\Sigma_{\rm d}$ (in $\S$4.2). We compare those distributions in star-forming regions against the rest of area in galaxies in $\S$4.3.

4.1 Distributions of Dust Mass and Luminosity as a Function of Illuminating Radiation Field

The dust mass distribution as a function of illuminating interstellar radiation field, $U$, strongly affects the integrated SED of a galaxy. Though we cannot measure the sub-resolution scale distribution (e.g., photo-dissociation region or Hii region structure), our maps allow us to measure the dust mass distribution as a function of $U$ on galactic scales. We measure and plot these distributions in the top panels of Figure 2 (for the Magellanic Clouds at 13 pc resolution) and the top panels of Figure 3 (for all targets at $167$ pc resolution).

Following Dale et al. (2001), we write these distributions as $dM_{\rm d}/d{\rm log}_{10}U$, where $dM_{\rm d}$ is the amount of dust mass heated by $U$ within a range of $\pm 0.5~{}d{\rm log}_{10}U$. We choose $d{\rm log}_{10}U=0.01$ dex, which is fine enough to capture the details of the distribution function, but not too fine that noise dominates the plots.

We normalize $dM_{\rm d}$ by the total dust mass, $M_{\rm Tot}$, in each galaxy, and define the mass fraction as $dm_{\rm d}=dM_{\rm d}/M_{\rm Tot}$. Then, the normalized distribution function is

[TABLE]

We interpret this normalized distribution function, $dm_{\rm d}/d{\rm log}_{10}U$, as the fraction of the total dust mass per dex of $U$, calculated in a bin with width of $d{\rm log}_{10}U$ and centered at $U$. Less formally, this is the normalized probability density of dust mass as a function of $U$.

Uncertainty: To assess the uncertainty in the measured distributions, we calculate the scatter in the distribution function across a set of 100 realizations of the fitting results (following Gordon et al., 2014; Chiang et al., 2018). We indicate this uncertainty by a shaded area in Figures 2 and 3, which shows the $16^{\rm th}{-}84^{\rm th}$ percentile across all realizations. The uncertainty appears larger at 167 pc resolution because lower number of pixels leads to more statistical noise.

To create the realizations, we use the relative likelihoods for each point in our $\Sigma_{d}{-}T_{d}$ grid space, and randomly draw values from that grid weighted by those relative likelihoods. We repeat this process 100 times for each pixel in each galaxy. This yields 100 realizations of the maps of $T_{d}$ and $\Sigma_{d}$ for each galaxy.

We construct the distribution functions as described above for each realized map. In each bin of the distribution function, we define the $\pm 1\sigma$ uncertainty as the difference between the $84^{\rm th}$ and $16^{\rm th}$ percentile value across all realizations in that bin.

Parameters of the distribution: For each distribution, we measure three quantities described as follows.

$U_{\rm peak}$, defined as the value of $U$ that corresponds to the peak of the distribution of mass or luminosity. This is the most common value in the distribution function. It often, but not always lies, near the center of the distribution. 2. 2.

$U_{\rm med}$, defined as the median value of $U$ in the distribution of mass or luminosity. 3. 3.

$\mathcal{W}_{U}$, defined as the logarithmic width in dex of $U$, that covers the 16th-to-84th percentile of mass or luminosity sorted by $U$.

We indicate $U_{\rm peak}$ and $U_{\rm med}$ for each galaxy by a circle and a tick above the distributions in Figures 2 and 3, and the width as the horizontal line, color coded by galaxy. We record their values in Table 3.

4.1.1 Distribution of Dust Mass as a function of Radiation Field

In practice, we calculate $dm_{d}/d{\rm log}_{10}U$ by summing the dust mass within each 0.01 dex-wide bin of log${}_{10}U$. Then, we divide the mass in each bin by the total dust mass in the galaxy and by the 0.01 dex of bin width. For bins of $U$ that have no mass in it, we assign the value of $dm_{d}/d{\rm log}_{10}U$ in those bins through linear interpolation. Before plotting, we apply Gaussian smoothing with width of five bins. We plot these distributions in the top left panels of Figures 2 and 3.

Because our estimate of $U$ depends directly on our best-fit of dust temperature, $T_{\rm d}$, via Equation 10, our measured distribution function can also be converted to the distribution of mass as a function of $T_{\rm d}$ via

[TABLE]

In the top left panels of Figures 2 and 3, we show these distribution functions in SMC and LMC at a common 13 pc resolution and in four targets at a common 167 pc resolution.

At 167 pc resolution, M31 has log${}_{10}U_{\rm peak}\approx-0.01$ dex (corresponds to $T_{\rm peak}\approx 18$ K). This is very close to the Solar neighborhood value of log${}_{10}U_{\odot}=0$ dex. But, as the top left plot in Figure 3 shows, $U_{\rm peak}$ in M31 lies at the upper end of a wide distribution (i.e. $U_{\rm peak}>U_{\rm median}\approx-0.14$). Therefore, M31 includes a large amount of mass illuminated by radiation fields below the Solar Neighborhood value. This agrees with the well-known result that the integrated SED in M31 indicates cooler dust temperatures than those in the other Local Group galaxies (e.g., Haas et al., 1998; Groves et al., 2012; Smith et al., 2012; Draine et al., 2014). As Figure 1 shows, much of this cooler material lies in the outer part of the galaxy.

At $167$ pc resolution, M33 shows the lowest log${}_{10}U_{\rm peak}$ in our sample222Keep in mind that we use a constant value of the dust emissivity index, $\beta=1.8$, for all galaxies. Since $T_{\rm d}$ and $\beta$ are anti-correlated, a lower value of $\beta$ in M33 ($\beta=1.5$; Xilouris et al., 2012) would lead to higher $U_{\rm peak}$. Conversely, a higher value of $\beta$ in M31 ($\beta\gtrsim 1.9$; Smith et al., 2012) would lead to a lower $U_{\rm peak}$. ($-0.52$ dex), but this peak lies towards the low end of a wide distribution (i.e. $U_{\rm peak}<U_{\rm median}\approx-0.33$). M33 has both more mass at high $U$ and more mass at low $U$ compared to M31 (shown as their horizontal lines), consistent with strong dust temperature gradient (e.g., Xilouris et al., 2012) and extended gas disk in that galaxy (e.g, Koch et al., 2018). Again, Figure 1 shows that the temperature and luminosity have a stronger concentration towards the inner part of the galaxy than the dust mass.

At $167$ pc resolution, both Magellanic Clouds show narrower distributions and higher $U_{\rm peak}$ than M31 and M33. The LMC has the highest log${}_{10}U_{\rm peak}$ (0.30 dex) in our sample, while the SMC has log${}_{10}U_{\rm peak}$ of 0.16 dex. These values are close to their log${}_{10}U_{\rm med}$ values (0.28 dex for LMC and 0.15 dex for SMC). Hence, both galaxies have most of their mass above the Solar neighborhood value ($U_{\odot}=1$). The LMC includes a significant “tail” toward high $U$, indicating a substantial fraction of its dust mass in high radiation field regions like that around 30 Doradus (see §4.3 and the maps in Figure 1). Compared to the two spirals (M31 and M33), the maps of the Magellanic Clouds show less mass in an extended, cool disk (Figure 1).

4.1.2 Distribution of Luminosity as a function of Radiation Field

We calculate the distribution of equilibrium IR luminosity as a function of $U$, $d\mathcal{J}_{\rm IR}/d{\rm log}_{10}U~{}(U)$, analogous to the calculation of mass distribution as a function of $U$. We normalize the distribution by the integrated luminosity for the whole galaxy, $\mathcal{J}_{\rm eq,Tot}$. The resulting distribution shows the fraction of the equilibrium luminosity per dex of log${}_{10}U$. We plot these distributions in the top right panels of Figures 2 and 3.

As one might expect, given the dependence of $\mathcal{J_{\rm IR}}$ on $U\propto T_{\rm d}^{5.8}$ (Equation 11), higher fraction of the luminosity distributions shift towards higher $U$ compared to the mass distributions. This strongly affects $U_{\rm peak}$ in M33, which has the lowest $U_{\rm peak}$ and $U_{\rm med}$ by mass but has the highest $U_{\rm peak}$ by luminosity (0.44 dex). The other galaxies show modest shifts in $U_{\rm peak}$ and $U_{\rm med}$, but the shape of their luminosity distributions around the peak change significantly (compared to that in the mass distribution). In both Magellanic Clouds, the luminosity distribution shows more prominent “bumpy” features at high $U$ than in the mass distribution. These reflect large contributions to the luminosity, but smaller contributions to the mass, from hot star forming regions (see §4.3).

At 167 pc resolution, M31 has the lowest log${}_{10}U_{\rm peak}$ (0.04 dex) and lowest log${}_{10}U_{\rm med}$ (0.03 dex), again, close to the Solar neighborhood value. Compared to the mass, less luminosity comes from the low temperature, extended part of the galaxy. This is evident where log${}_{10}U_{\rm med}$ for luminosity is $0.03$, but only $-0.14$ for the mass distribution (i.e. lower value of $U_{\rm med}$ means low-$U$ regions contribute more toward the overall distribution). The prominent low-$U$ feature in the mass distribution appears dramatically suppressed in the luminosity distribution. Meanwhile, the central part of the distribution ($U\sim 1$), that maps to the star-forming 10 kpc ring in M31 (Figure 1), increases in prominence. The high-$U$ bulge region contributes much more luminosity than mass, but still represents only a small fraction of the galaxy’s total luminosity (1.6%; see §4.3 for detailed discussions).

M33 shows the most striking change between the mass and luminosity distributions. Its $U_{\rm peak}$ shifts from the lowest by mass ($-0.52$ dex) to the highest by luminosity ($0.44$ dex). Its $U_{\rm med}$ also increases by 0.60 dex. This reflects the strong dust temperature gradient (Xilouris et al., 2012), which creates a wide range of $U$ and $T_{\rm d}$, along with the extended gas and dust surface density distributions (Koch et al., 2018). As the maps in Figure 1 show, the mass in M33 extends out to large radii, while the luminosity appears more centrally concentrated, coincides with active star formation (§4.3).

We note that in SMC and LMC, there is a fraction of IR luminosity that comes from dust with $T_{d}>30$ K in the star forming regions (high end tail in the top right panel of Figure 2). Even though warm dust component can contribute towards these hot regions, it only contributes a negligible fraction in the mass distribution and lies outside 68% of the mass distribution (top left panel of Figure 2). Hence, the accuracy of our modeling for this warm dust component does not affect the main results of this paper.

4.1.3 The Width of the Distributions

At 167 pc resolution in the Magellanic Clouds and M31, $68\%$ of both the mass and IR luminosity are concentrated within a $\lesssim 0.6$ dex range ($\pm 0.3$ dex, or about a factor of two). In the SMC, the distribution of luminosity appears even narrower, with $68\%$ coming from a $0.37{-}0.43$ dex wide range.

M33 yields the widest distributions of both mass (1.04 dex) and luminosity (0.82 dex) as a function of $U$. As discussed above, these wide distributions reflect the structure of M33 that is visible in Figure 1. High $\Sigma_{\rm d}$ extend into the outer disk, while $T_{\rm d}$ tends to drop with radius. Meanwhile, in the inner part of the galaxy, multiple knots of star formation activity exhibit both high $T_{\rm d}$ and $\Sigma_{d}$.

At $13$ pc resolution, available only for the Magellanic Clouds (Figure 2), the distributions appear wider than at $167$ pc resolution. These distributions at 13 pc resolution show more material at low radiation fields. This could be expected if the convolution to coarser resolution blurs together warm and cool dust. At their coarser resolution, the light from nearby warmer dust would make it harder to isolate the cool component, leading to the lower amount of low-$U$, low-$T_{\rm d}$ materials at $167$ pc resolution. This idea of relatively “hidden” cool dust, masked by the presence of more luminous material nearby, has appeared many times in the literature (e.g., Galliano et al., 2005). We explore this effect quantitatively in $\S$6.1.

4.2 Distributions of Dust Mass and Luminosity as a Function of Dust Mass Surface Density

We also build the distributions of mass and luminosity as a function of the dust mass surface density, $\Sigma_{d}$. We plot these $dm_{d}/d{\rm log}_{10}\Sigma_{d}$ and $d\mathcal{J}_{\rm IR}/d{\rm log}_{10}\Sigma_{d}$ in the lower panels of Figures 2 and 3. We construct these distributions in the same way as those treating $U$ as the independent variable (in $\S$4.1), but here, we bin it by log${}_{10}\Sigma_{d}$ instead of log${}_{10}U$.

As in $\S$4.1, we characterize the distributions with three parameters, $\Sigma_{d}$ at the peak of the distribution ($\Sigma_{\rm peak}$), the median of the distribution ($\Sigma_{\rm med}$), and the logarithmic width of the distribution, $\mathcal{W_{\Sigma}}$, defined as the 16th-to-84th percentile range of the distribution. We show these $\Sigma_{\rm peak}$ and $\Sigma_{\rm med}$ as the dots and ticks above the distributions in Figures 2 and 3.

4.2.1 Mass Distribution as a function of Dust Mass Surface Density

The bottom left panel in Figure 3 shows the dust mass distribution for our targets at 167 pc resolution. Dust mixed with atomic gas, molecular gas, and even ionized gas are all contribute to these histograms, weighted by the local gas-to-dust ratio. In that sense, this plot also shows the column density distribution of the whole ISM across the Local Group.

Almost all dust mass in all four targets lies in the range $\Sigma_{d}\sim 0.01{-}0.1$ $M_{\odot}$ pc*-2*. For a Galactic gas-to-dust ratio (GDR) of $\sim 150$, this range equates to a gas mass surface density of $\sim 1{-}15~{}M_{\odot}$ pc*-2*. The metallicities of the LMC, SMC, and M33 are all lower than that in the Milky Way (e.g., Russell & Dopita, 1992; Rosolowsky & Simon, 2008), so the associated range of gas mass surface densities should in fact be more like $\sim 3{-}50~{}M_{\odot}$ pc*-2*. This range is consistent with the fact that the ISM in all of these galaxies is known to be dominated by atomic gas across the disk, with a few dense regions (e.g., the inner part of M33 and the LMC’s ridge region) locally dominated by molecular gas (e.g., Druard et al., 2014; Wong et al., 2011).

The LMC and M33 show similar mass distributions. They exhibit the highest $\Sigma_{\rm peak}$ of $\approx 0.045~{}M_{\odot}$ pc*-2* in our sample (their $\Sigma_{\rm peak}\approx\Sigma_{\rm med}$) and widths of $0.5{-}0.7$ dex. These two dwarf spirals have comparable metallicity, $Z\approx 0.5~{}Z_{\odot}$ (Russell & Dopita, 1992; Rosolowsky & Simon, 2008), and the distribution appears consistent with dust being mostly mixed with high column density of Hi in both targets.

The SMC shows low $\Sigma_{\rm peak}$ and low $\Sigma_{\rm med}$ of $0.02$ and $0.03$ $M_{\odot}$ pc*-2*, respectively. The SMC has a lower metallicity and higher GDR of $\sim 1200$, compared to those in M33 and the LMC ($\approx 380$; Roman-Duval et al., 2014). The lower $\Sigma_{\rm peak}$ in the SMC could thus be expected if the gas mass surface density in SMC were comparable to that in LMC. In reality, the SMC’s elongated structure along the line of sight (Stanimirovic et al., 1999; Staveley-Smith et al., 2003) leads to higher Hi column densities than one finds in the other Local Group targets. These effects combine to produce the multi-components distribution (bumps) that mostly overlaps the other targets in the bottom left panel of Figure 3. We emphasize that this overlap is at least partially a coincidence due to the SMC’s orientation.

M31 also shows low $\Sigma_{\rm peak}$ and $\Sigma_{\rm med}\sim 0.02$ $M_{\odot}$ pc*-2*, again indicating a large amount of mass in an extended disk with relatively low column densities and a GDR that increases with radius (e.g., Draine et al., 2014). Note that M31’s appearance in the plot depends on the substantial correction that we have applied to account for M31’s high inclination ($77^{\circ}.7$; Corbelli et al., 2010).

4.2.2 Luminosity Distribution as a function of Dust Mass Surface Density

The distribution of luminosity, $\mathcal{J_{\rm IR}}$, shifts toward higher $\Sigma_{d}$ compared to the distribution of mass. The contribution of mass from low $\Sigma_{d}$ is suppressed in the luminosity distribution, changing the shape of distribution in M31 and M33 from symmetric, flat/roundish top, to have narrow peak with tail at low end (positively skewed). The main peak in the luminosity distribution for M31 mostly comes from the star-forming ring and less from the low $T_{d}$ in the outer disk. The tail towards low $\Sigma_{d}$ shows the influence of the hot bulge.

In the LMC, the distinct components (bumps) in the mass distribution become more pronounce in the luminosity distribution, reflecting the fact that warm dust in high-$\Sigma_{d}$ star-forming regions contribute a lot to the IR luminosity. In the SMC, the distribution changes from negatively skewed in mass to a more symmetric in luminosity, because the high mass surface density regions also tend to have warmer temperature that emits a large fraction of the luminosity.

Because IR luminosity proportional to $\Sigma_{d}T_{d}^{5.8}$, the difference between the mass and luminosity distributions implies that temperature correlates with $\Sigma_{d}$, i.e. part of the luminosity distribution from low $\Sigma_{d}$ is suppressed (because they also have lower $T_{\rm d}$), while luminosity distribution from high $\Sigma_{d}$ is enhanced (because they also have higher $T_{\rm d}$). In $\S$5.1 we show that this appears to be the general case in high signal-to-noise regions, though with some notable exceptions.

4.2.3 Widths Comparisons

At 167 pc resolution, the SMC and LMC have similar width in the distributions of mass and luminosity ($\sim 0.7$ dex), and only slightly larger than the width for M31 and M33 ($\sim 0.5{-}0.6$ dex). This means the SMC and LMC have the largest variations in $\Sigma_{d}$ where 68% of the mass is distributed. As in the mass distribution, the width of luminosity distribution at 167 pc is smaller than that at 13 pc resolution (Table 3) because of convolution between regions with low and high mass surface density.

The shapes of the mass distributions somewhat resemble the distributions of molecular gas mass as a function of $\Sigma_{d}$ in nearby galaxies seen at similar spatial resolution by Sun et al. (2018). They also found signature of multiple components, usually related to distinct dynamical regions (bulge vs. disk). The distributions that we find show more width and more components, perhaps reflecting the wider mixture of environments (e.g. spiral arms and star-forming regions) that is incorporated into our distribution.

In §4.3, we will see that the distributions for any isolated region may appear roughly lognormal in shape. This would agree with the observation that the column density distribution of diffuse neutral gas density in the Milky Way and other Local Group galaxies has a lognormal shape (Berkhuijsen & Fletcher, 2008, 2015; Corbelli et al., 2018). The same shape also emerges from simulations (Wada & Norman, 2007). Making a more rigorous measurement of the shape of the distribution after controlling for regional variations is a subject for future work.

4.2.4 Implications

These mass and luminosity distributions as a function of $\Sigma_{d}$ have two implications: low infrared optical depth and no evidence for opaque Hi, as described below.

Low Infrared Optical Depth: The dust mass surface densities in our targets almost never exceed $\Sigma_{\rm d}\sim 1$ M⊙ pc*-2*. Given our adopted value of $\kappa$ at 160 µm, this implies that at our resolutions, the optical depth throughout our targets at 100 and 160 µm will be $\lesssim 1\times 10^{-2}$ and $5\times 10^{-3}$, respectively. This validates our assumption of optically thin dust that is used in the model (i.e. approximating Equation 1 by Equation 2).

The opacity at $100{-}160$ µm appears small at both $13$ pc (for the Magellanic Clouds) and $167$ pc (for all targets). The importance of pressure from reprocessed IR emission scales as the infrared opacity, $\tau_{IR}$ (e.g., see Thompson et al., 2005), and should exceed the UV-driven radiation pressure when $\tau_{\rm IR}>1$. Thus, we expect radiation pressure from long wavelength infrared photons will contribute an insignificant amount of feedback over essentially all of the area in these targets.

The highest $\Sigma_{d}$ that we observe in the LMC (Figure 2) imply significant opacity at shorter IR wavelengths (e.g., $\lambda\sim 20$ µm). We do not consider emission at these wavelength. However, given the surface densities that we observe and the fact that star forming regions often produce significant emission at these mid-IR wavelengths, our observations are consistent with a significant effect from reprocessed infrared photons in the hottest, highest density parts of our targets (Lopez et al., 2011, 2014). To see this, note that at $13$ pc resolution, the LMC by luminosity includes a small, but noticeable, contribution at $\Sigma_{d}\sim 1{-}10$ M⊙ pc*-2* (bottom right panel of Figure 2). This $\Sigma_{d}\sim 1$ M⊙ pc*-2* implies $\tau_{\rm 160\micron}$ of only $\sim 4\times 10^{-3}$. If we extrapolate to $\sim 20$ µm, though, $\tau_{\rm 20\micron}\sim 0.15{-}0.3$ (the extrapolation using $\beta\sim 1.8{-}2$ is quite approximate but not unreasonable, see Draine 2011). Then, even at the extreme tail of the LMC distribution, we should see a non-trivial, and perhaps even dominant, contribution to the radiation pressure from IR-reprocessed photons.

No Evidence for Hidden Opaque Hi: The modest $\Sigma_{\rm d}$ values that we observe and the smooth appearance of our $\Sigma_{\rm d}$ maps provide evidence against the existence of the patchy, high column density, opaque Hi clouds posited by Braun et al. (2009) and Braun (2012). Braun presented opacity-corrected Hi column density maps in three of our targets; the LMC, M31, and M33. The 21-cm maps that were used to create these Hi maps have comparable angular resolution to our $\Sigma_{\rm d}$ maps. In Braun’s maps, the high gas column density features ($N_{\rm H}\sim 10^{23}$ cm*-2*) appear as small patches sprinkled throughout the maps.

In contrast, our $\Sigma_{\rm d}$ maps do not have a similar patchy appearance, nor do the values of $\Sigma_{\rm d}$ that we find support the presence of gas with column densities approaching $10^{23}$ cm*-2*, scattered through the Local Group at this resolution. For the LMC, M33, and M31, we expect the gas to dust ratios of $\sim 100{-}500$ (Leroy et al., 2011; Roman-Duval et al., 2014). In those cases, the $\Sigma_{\rm d}$ that we find corresponds to column densities of $<5\times 10^{22}$ cm*-2*. In order for the opaque Hi features posited by Braun (2012) to exist but not appear in our $\Sigma_{\rm d}$ maps, they would need to be unusually dust poor (i.e. gas to dust ratio of $>1,000$). To our knowledge, there is no plausible mechanism to produce dense, opaque atomic gas clouds that are preferentially depleted in dust.

Our results do not rule out 21-cm line opacity playing an important role in galaxies, even in these galaxies. More modest effects that operate smoothly across the maps would be hard to distinguish from variations of the gas-to-dust ratio or the presence of CO-dark molecular gas. This issue certainly remains an open topic. Our maps simply do not support with specific structure found by the Braun’s analysis.

4.3 Regional Comparisons

As we see in Figures 2 and 3, there are clear differences in the shape of distribution between mass (left panels) and luminosity (right panels). There are also multiple components in some of the distributions. Comparing the distributions to the maps in Figure 1 suggests that some of these features arise from prominent star-forming regions. Features associated with ongoing massive star formation appear particularly prominent in the luminosity maps of all four targets.

To make a quantitative connection between the structure in the maps and the distribution functions, we divide each galaxy into “star-forming” and “non star-forming” regions and build separate distribution functions for each. We normalize these mass and luminosity distributions by the total mass and total luminosity for the whole galaxy (as in §4.1 and §4.2). For M31, we also build separate distributions for the hot, low $\Sigma_{d}$ central region (inside 1 kpc from the nucleus).

Selection of Regions: We define “star-forming” (SF) regions as areas with Spitzer MIPS (Werner et al., 2004; Rieke et al., 2004) 24 µm intensity $>1$ MJy sr*-1* (shown as the blue regions in the right panels of Figures 4 and 5). This intensity is equivalent to $\approx 7\times 10^{-4}~{}M_{\odot}$ yr*-1* kpc*-2* (using the calibration from Calzetti et al., 2007). We choose this particular value because it selects the bright regions in the star-forming ring of M31 (Tabatabaei & Berkhuijsen, 2010) and star-forming complexes in the Magellanic clouds and M33 (Relaño & Kennicutt, 2009).

This cutoff also selects the bulge of M31, which shows hot dust (in Figure 1) that Groves et al. (2012) and Smith et al. (2012) argued to be heated by the older stellar population. Therefore, we build separate distributions for this central region (defined as the area within 1 kpc from the nucleus, based on a position angle of $37^{\circ}.7$ from Corbelli et al., 2010).

Resulting Distributions: We show the mass and luminosity distribution functions for star-forming regions (blue) and the rest of area in the galaxy (red) in Figures 4 and 5 (with a log-scale in the vertical axis).

In all four galaxies, both the mass and the luminosity distributions from star-forming regions appear distinct compared to the rest of area in the galaxy. The star forming regions tend to show narrower distributions with a few prominent peaks. As expected, in all four galaxies, the star-forming regions tend to have warmer dust temperature than the rest of area in the galaxy.

Those figures also show that multi-component distributions in Figures 2 and 3 are a result from superimposing physically distinct regions. Star-forming regions, which also tend to be the regions rich in molecular gas, contribute many of the high-$\Sigma_{d}$ and high-$U$ features in the full-galaxy distributions. Because the star forming regions have both high $T_{d}$ and high $\Sigma_{d}$, they appear even more prominent in the luminosity distribution.

To be concrete, for our adopted definition, the fraction of dust mass that resides in the star forming regions is 10.0% in the SMC, 28.4% in the LMC, 10.2% in M31, and 13.9% in M33. The fraction of luminosity from the star-forming regions is higher, i.e. 24.3% in SMC, 51.4% in LMC, 19.4% in M31, and 44.8% in M33.

The bulge in M31 has higher dust temperature than the star-forming ring or the rest of the galaxy ($\approx 30$ K vs. 17 K), but it only captures very small fraction of the total dust mass ($\approx 0.08\%$). However, this central region (orange histogram in Figure 5) becomes more prominent in the luminosity distribution, with luminosity fraction of $1.6\%$. This more prominent contribution in the luminosity distribution, in general, is also true for star-forming regions.

5 Resolved Correlations Among Dust Parameters

5.1 Correlation between Dust Temperature and Mass Surface Density

Both the maps (Figure 1) and the distributions (Figures 2 to 5) suggest a correlation between $\Sigma_{d}$ and $T_{d}$, so that high mass surface density regions also tend to be hotter and more luminous. Figure 6 shows this relationship directly, where we plot dust temperature against the dust mass surface density pixel-by-pixel at 13 pc resolution for the Magellanic Clouds and 167 pc resolution for M31 and M33.

At low signal to noise, we expect the correlated uncertainties to drive an apparent anti-correlation between the best fit $\Sigma_{d}$ and $T_{d}$. This reflects that for higher temperature, less mass will be required to produce any given intensity (Equation 5). To suppress this effect, we plot pixels with signal-to-noise less than $3$ in either $T_{\rm d}$ or $\Sigma_{d}$ as gray points and contours. For the rest of this section, we focus our attention on the blue points and contours, which have $S/N>3$ in both quantities and should be less affected by noise.

Because the maps contain many individual pixels, we use the Gaussian kernel density estimate in Scipy (with bandwidth estimate following Scott, 1992) to calculate the probability density of data points in $200\times 200$ grids between the minimum and maximum values of those data points. The contours in Figure 6 mark the probability density of 0.005, 0.01, 0.05, and 0.1 data points per grid.

For the SMC, LMC, and M33, the red points show the median temperature in 0.1 dex bin of log${}_{10}\Sigma_{d}$, in regions where log${}_{10}\Sigma_{d}>-1.5$. Even though the Spearman (1904) rank correlation coefficient, $r_{s}$, is low for all blue points (S/N $>3$), the weak correlation between temperature and dust mass surface density is highly significant because the $p$-value is almost zero. This means warm regions also tend to have high mass surface densities. As we saw in the last section, star forming regions in the maps tend to have both high mass surface density and high temperature.

At some degrees, a correlation between $T_{\rm d}$ and $\Sigma_{d}$ should be expected on large scales for dust internally heated by star formation. Following Schmidt (1959) and Kennicutt (1998), gas with high mass surface density forms stars more rapidly. If the light from young stars is reprocessed into equilibrium IR emission, then the emergent luminosity should track this star formation. If we approximate $\Sigma_{\rm SFR}\propto\mathcal{J_{\rm IR}}\equiv\Sigma_{d}T_{d}^{5.8}$ and $\Sigma_{\rm gas}\propto\Sigma_{d}$, then the correlation between $\Sigma_{d}$ and $T_{d}^{5.8}$ appears when the $\Sigma_{\rm SFR}-\Sigma_{\rm gas}$ relationship is steeper than linear. Indeed, these galaxies are all mostly Hi-dominated, and the $\Sigma_{\rm SFR}-\Sigma_{\rm HI}$ scaling relation does have a super-linear slope (e.g., Bigiel et al., 2008; Schruba et al., 2011), but also has an important dependency on other parameters (Leroy et al., 2008). This one-to-one correlation between $\Sigma_{d}$ and $T_{d}^{5.8}$ is shown as the black line in Figure 6 (with an arbitrary normalization).

This $\Sigma_{d}\propto T_{d}^{5.8}$ scaling would be expected for dust heated by star formation. But our sample also includes several cases where this should not be a good approximation. As mentioned several times, the bulge of M31 has hot, low surface density dust. This stands out in the bottom left panel of Figure 6. This dust has been shown to be externally heated by the old stellar population (Groves et al., 2012; Draine et al., 2014).

At $13$ pc resolution, the approximation of internal heating may break down due to time evolution of star forming regions. Following Schruba et al. (2010) for M33, at such high resolution the heating sources and gas peaks resolve into discrete features. This effect can also be seen in the Magellanic Clouds (Jameson et al., 2016). In this case, the local $\Sigma_{d}$ will no longer be directly correlated to the amount of heating.

5.2 Intensity at 500 µm and Dust Mass Surface Density

With powerful sub-millimeter telescopes like ALMA, it has become common to estimate the dust mass from only one or a few measurements on the Rayleigh-Jeans tail of the SED. This approach assumes a dust temperature, either implicitly or explicitly. Here, we compare $\Sigma_{d}$ from our five bands fitting, which leverages information from the $100{-}500$ µm SED to the monochromatic intensity at $500$ µm. This offers a check on the reliability of sub-mm intensities as a dust mass tracer.

Figure 7 shows the correlation between resolved $I_{500\micron}$ and $\Sigma_{d}$ over three orders of magnitude in both axes. We only select pixels above $2\sigma$ uncertainty for both $I_{500\micron}$ and $\Sigma_{d}$. Contours show the density of data points, generated using the Gaussian kernel density estimate with a kernel size of $\approx 0.02$ dex. In each panel, the diagonal lines show relations with power law index of 1 with normalization set by the median of $\Sigma_{d}/I_{500\micron}$ measured for each galaxy. The value of this median ratio is $\Sigma_{d}/I_{500\micron}\approx$ $1.91_{-0.43}^{+0.51}\times 10^{-2}$ for the LMC, $2.05_{-0.53}^{+0.71}\times 10^{-2}$ for the SMC, $2.84_{-0.65}^{+0.94}$ for M31, and $2.70_{-0.87}^{+1.31}$ for M33, where $\Sigma_{d}$ is in units of $M_{\odot}$ pc*-2* and $I_{500\micron}$ is in units of MJy sr*-1*. The lower and upper limits are the 16th and 84th percentiles, respectively. Among all of our targets, the variation of this ratio roughly agrees with the simple expectation based on our measured variations in $T_{d}$.

Both contours of M31 and M33 show a shallow slope in Figure 7 (compared to one-to-one relation). These reflect the temperature gradients discussed above (Xilouris et al., 2012). In these galaxies, taking $I_{500\micron}$ to linearly trace $\Sigma_{d}$ would give a bias so that the dust appears to be more centrally concentrated than it actually is. In other words, the temperature gradient causes the inner part of M33, and the star-forming ring and bulge of M31, to glow more brightly than the rest of the galaxy, while the cooler outer regions are fainter. This effect appears less pronounced in the Magellanic Clouds at high resolution.

Overall, this exercise confirms that $I_{500\micron}$ predicts $\Sigma_{d}$ within an accuracy of about $50\%$ for our sample, but with systematic biases in the two disk galaxies with temperature gradients. We emphasize that this only tests the ability to recover the dust surface density. We have taken no account here for variations in the gas-to-dust ratio.

5.3 Infrared Color and Temperature

Finally, we compare the 100-to-250 µm or 100-to-350 µm infrared color to the $T_{d}$ from our five band fits. This allows us to check how well our model reproduces this IR color and to verify how well we would have predicted $T_{d}$ from only a single band ratio. We focus on $I_{100\micron}/I_{250\micron}$ and $I_{100\micron}/I_{350\micron}$ as colors because those wavelengths bracket the peak of the dust SED at $\sim 160$µm.

Figure 8 shows $T_{d}$ as a function of IR color for data with $S/N>10$ in both color and temperature. Contours indicate data density as in §5.2. We also plot the relation between temperature and colors as expected from the modified blackbody models with $\beta=$ 1.4, 1.8, and 2.2.

The data broadly follow the model, although with some notable deviations. The plots suggest that LMC and M33 may be better described by $\beta$ closer to $\sim 2$ than our adopted $\beta=1.8$. Both the SMC and the LMC show points covering a range of $\beta$, with the regions with redder IR color are more consistent with higher $\beta$, while the rest of the pixels is in good agreement with our adopted value of $\beta$. Note that despite the shift in $\beta$, these features are as likely come from multiple dust populations unresolved by our observation (mixing of multiple populations within the beam) and the influence of out-of-equilibrium heating at short wavelengths, so that they appear as variations in dust properties (changing $\beta$). These reflect deviations from our model. The strength of this feature in the SMC appears consistent with previous observations that hot, small grains contribute heavily to the SED in this galaxy (e.g., Bot et al., 2004).

Lastly, we note that in the range of $30<T_{d}<40$ K, the modified blackbody models still show a curvature in the color${-}$temperature relation (black curves in Figure 8). This means we can still constrain the dust temperature well within that temperature range. A vertical slope would mean the temperature loses its sensitivity to the IR color (i.e. approaching Rayleigh-Jeans tail), and this is clearly not the case in our study. Therefore, the high-temperature tail in the distributions (Figure 2) is real.

6 The Dependence of Dust Mass and Temperature on Physical Resolutions

6.1 Dust Temperature as a Function of Physical Scale

We study four of the closest star-forming galaxies at the diffraction limit of the best infrared telescope to date. This high physical resolution is not available for distant galaxies. Still, as discussed in §1, there remains considerable interest in using dust as a tracer of gas in more distant systems. Therefore, we explore the effects of physical resolution on the derived dust temperature and mass.

We choose the weighting methodology (described below) to create coarser resolution maps. This is different than convolving the SED first and then fitting that convolved SED. The reason behind this choice is we lose the background area gradually in the coarser resolution maps. This would make the covariance matrix calculation (Equation 7) becomes difficult and leading to non-robustness of the fitting results of the convolved SED (see Appendix A).

Method: Observing at coarse resolution tends to emphasize the sources of luminosity, rather than mass. To investigate the effects of resolution, we compare our high resolution results to what we would derive from luminosity (rather than mass) weighted averages at some coarser resolutions. We compare any results at coarser resolutions to the mass-weighted average at the finest resolution, because we consider the latter as a quantity closest to the true temperature from the majority of dust mass.

Specifically, we calculate the luminosity-weighted temperature, $T_{d,L}$, and the luminosity-weighted mass surface density, $\Sigma_{d,L}$, at resolutions coarser than our original maps. To calculate $T_{d,L}$, we begin with our highest resolution $T_{\rm d}$ map, weight that map by luminosity, convolve that map to lower resolution, and then, divide the result by the convolved luminosity map. Mathematically, this is expressed as

[TABLE]

where $\mathcal{J_{\rm IR}}$ is a quantity proportional to the equilibrium dust luminosity per unit area in the highest resolution map (as defined by Equation 11), $\ast$ denotes convolution, and $\mathcal{G}$ is the Gaussian kernel to convolve from the finest resolution to the target coarser resolution.

Similarly, we calculate the luminosity-weighted mass surface density, $\Sigma_{d,L}$, via

[TABLE]

In Equations 14 and 15, we calculate $T_{d,L}$ and $\Sigma_{d,L}$ for each pixel at coarser resolution. To represent the temperature for a single galaxy, we also calculate the mean value of temperature (weighted by mass) at each resolution via

[TABLE]

Similarly, the luminosity-weighted mean temperature for a single galaxy at coarser resolution is

[TABLE]

Here, the sum $\sum_{i}$ runs over all pixels inside the region of interest.

As we have seen in $\S$4.3, weighting by luminosity leads to higher $\langle T_{\rm d}\rangle$ because more weight is attached to hot, high luminosity star-forming regions. Because we use the luminosity-weighting to convolve the data to coarser resolutions, information below the resolution is irrevocably “washed out.” Comparing mass-weighted results at different resolutions shows how much this “sub-resolution” luminosity-weighting changes our inferences about the dust temperature and mass. This gives us an empirical estimate of how much the implicit luminosity-weighting in low resolution observations will wash out information on the true dust mass.

Results for temperature: In Figure 9, we plot $\langle T_{\rm d}\rangle_{M}$ (black dots) for each galaxy as a function of physical resolution. We also mark $\langle T_{\rm d}\rangle_{M}$ at the finest resolution as blue lines. This is the quantity that is closest to the true mass-weighted mean of temperature in a galaxy. Any other measurements should be compared to this line. Due to different physical sizes among galaxies, LMC and SMC are convolved to $7$ resolutions from $13$ pc to $500$ pc, while M31 and M33 are convolved to $5$ resolutions from $167$ pc to $2$ kpc.

In general, $\langle T_{\rm d}\rangle_{M}$ rises as the resolution becomes coarser. This has the expected sense, because convolving cool, low surface density regions with warm, high surface density regions would bias $\langle T_{\rm d}\rangle_{M}$ towards higher temperatures. We already saw an indication of this effect when comparing the distributions of mass by radiation field at $13$ pc and $167$ pc resolution for the Magellanic Clouds (Figures 2 and 3).

This trend of $\langle T_{\rm d}\rangle_{M}$ with spatial resolution is weakest in M31; the inferred mean temperature only changes by $<1$ K as we blur the maps by a factor of $\sim 10$ in resolution. We identify two likely causes. First, we see a weak correlation between temperature and mass surface density in Figure 6. Second, much of the variation in temperature occurs radially on kpc scales or larger, so that blurring together the distinct physical conditions is less of an issue in M31.

The SMC and M33 show the largest increases in $\langle T_{\rm d}\rangle_{M}$ with coarser spatial resolutions. The implied temperature in both galaxies increases by $\sim 2$ K as the resolution degrades by a factor of $\sim 10$. The LMC shows an intermediate case. We again attribute the magnitude of this resolution effect to the spatial structure and dynamic range in temperature variations in the galaxy. In general, the smaller, more star-formation dominated systems show more systematic effects, and more information loss, from blurring the data.

The red and blue horizontal lines in Figure 9 indicate the luminosity weighted temperature, $\langle T_{\rm d}\rangle_{L}$, and mass-weighted temperature, $\langle T_{\rm d}\rangle_{M}$, respectively, at the finest resolution. By construction, the black point lies on the blue line at the finest resolution. As we blur the galaxy towards a single resolution element, the black points should approach the red line. The difference between the red and blue lines shows the maximum effect of interchanging luminosity and mass weighting. This difference is $\sim 2{-}3$ K, or about $10\%$ change in $T_{d}$, or a factor of $\sim 2$ change in $U$.

Implication for the total mass measurement: Luminosity weighting tends to bias towards higher temperature. An overestimate of temperature will lead to an underestimate of the mass. To see this, consider a point on the Rayleigh-Jeans tail, for which $I_{\nu}\propto T_{d}~{}\Sigma_{d}$. As the luminosity-weighting biases toward higher temperature, it will also lead to an underestimate of the mass.

To first order, the mass will be biased low by the same fraction as the temperature is biased high. This is $\sim 10{-}20\%$ change over the range of resolutions that we study. However, this does not capture the effects of any distributions of mass or temperature beneath our finest resolution (13 pc for Magellanic Clouds and 167 pc for M31 and M33).

For the LMC, Galliano et al. (2011) considered exactly this effect and estimated that the mass inferred from the integrated SED is biased low by $50\%$. They included shorter wavelength bands than what we study here, which likely accounts for that stronger bias than what we found. The effect of luminosity weighting should be stronger for these shorter wavelength bands because they are more sensitive to the dust temperature.

In general, following Dale et al. (2001), accounting for distributions of temperatures in modeling the SED is crucial for a correct mass estimation. Emphasizing longer wavelengths, as we do here, somewhat diminishes this resolution effect. Even at 13 pc resolution, mass estimations in SMC and LMC are still less accurate because a significant portion of blending populations (unresolved by our study) must still occur even below that resolution.

6.2 Fits to the Integrated SED and Comparison with Resolved Mapping

We also fit the integrated SED, equivalent to unresolved observation. Even in the local universe, this is often the only type of measurement available for a galaxy. In particular, poor resolution options remain the main way to study high redshift galaxies (e.g., Magdis et al., 2010, 2012; Genzel et al., 2015) and there has been some controversy about how to best estimate dust mass in these galaxies. Here, following Galliano et al. (2011), we check how the temperature and dust mass derived from the integrated SED compared to those derived from the highly resolved maps.

We calculate the integrated SED by taking the mean surface brightness within the region of interest. We report the mean intensity for each target at each band in Table 4. We also note the averaging area, in steradians. To calculate the integrated SED in units of specific flux, we multiply the area by the intensities. Then, we fit the SED with two approaches; single-temperature and multi-temperature modified blackbodies. Figure 10 shows these integrated SEDs as red points, along with the best-fit models (curves). Uncertainty in the best-fit models, visualized by drawing 100 realizations in proportion to their likelihood, are shown as shaded areas.

In Table 4, we compare the dust temperature and total mass between galaxies from three different approaches; (1) integrated SED, fitted using single-temperature modified blackbody, (2) integrated SED, fitted using multi-temperature modified blackbody, and (3) the mass-weighted mean of temperature, $\langle T_{\rm d}\rangle_{M}$, calculated from the finest resolution map.

Single-Temperature Model: First, we use the same single-temperature modified blackbody model as adopted for the resolved observations. The best-fit models are shown as the solid green curves in Figure 10. Despite the simplicity of our model, it offers a good fit to the integrated SED of each galaxy. Mathematically, this is different than the luminosity-weighted mean of temperature, $\langle T_{\rm d}\rangle_{L}$, because in the integrated SED, we take average of intensity without any weighting. Therefore, we do not expect $\langle T_{\rm d}\rangle_{L}$ to be the same as the temperature derived from the integrated SED fitting (red vs. green lines in Figure 9).

The expected value of $T_{\rm d}$ from fitting these integrated fluxes are shown as the green lines in Figure 9. Dust temperature from fitting the integrated SED is higher than the mass-weighted temperature of the most highly resolved observations. This means the unresolved observations would overestimate the temperature of the majority of the dust by $\lesssim 3$ K or a factor of $\lesssim 17\%$ for 18 K temperature. If using a single flux at 500 µm to measure the dust mass (Equation 1 in Scoville et al., 2014), this translates to an underestimate of dust mass also by $\lesssim 17\%$ for an object with 18 K dust temperature.

Multi-temperature Model: We also fit the integrated SED by adopting a sub-resolution multi-temperature distribution in the form of (for $\alpha\neq 1$; Draine & Li, 2007)

[TABLE]

where $\delta$ is Dirac’s delta function. Here, a fraction $\gamma$ of the dust mass within a resolution element is heated by a power law distribution field $U^{-\alpha}$ between $U_{\rm min}$ and $U_{\rm max}$, while the rest is heated by the ambient ISRF of $U_{\rm min}$.

We use the same grid-based fitting procedure as before (§3.2). The grid for $\alpha$ is from $1.1$ to $5.1$ with an increment of $0.2$, the grid for log${}_{10}\gamma$ is from $-3$ to [math] with an increment of $0.2$, and the grid for log${}_{10}U_{\rm min}$ from $-1.5$ to $1.5$ with an increment of $0.2$. The grids for dust temperature and mass surface density are the same as in the single temperature modified blackbody model (§3.2). The value of $U_{\rm max}$ is fixed at $10^{3}$. The adoption of this $U_{\rm max}$ value does not affect the outcome of fitting parameters.

The best-fit models are shown as the dashed purple curves in Figure 10. The expected values for $U_{\rm min}$, $\alpha$, and $\gamma$ are listed in Table 5. The value of $U_{\rm min}$ represents the ambient ISRF that illuminate the majority ($1-\gamma$) of dust mass. Hence, it is not surprising that it is highest in the LMC, in general agreement with the distribution of dust mass as a function of $U$ ($\S$4.1). The rest of dust mass (a fraction of $\gamma$) is heated by ISRF in the range between $U_{\rm min}$ and $U_{\rm max}$.

The value of $\alpha$ shows the relative contribution of dust mass at the high-end tail of $U$, compared to that in the low-end tail of $U$, where $U_{\rm min}<U<U_{\rm max}$. A shallower slope (lower $\alpha$), such as in the SMC and M33, means the low-end tail of $U$ (originated from star-forming regions) contributed more to the dust mass spectrum between $U_{\rm min}$ and $U_{\rm max}$. Hence, this multi-temperature fit is not only adding free parameters, but also physically more reasonable than a single-temperature fit of modified blackbody.

Following Draine & Li (2007), we calculate the mass-weighted mean of $U$ via

[TABLE]

Then, we convert $\langle U\rangle_{M}$ to the dust temperature by inverting Equation 10.

Comparison between methods: In SMC and M33, the temperature derived from integrated SED multi-temperature is lower than those derived from the integrated SED single-temperature model by $\sim 1-2$ K, while the opposite happens in M31. The best-fit temperature from two methods above matches in LMC. Compared to $\langle T_{\rm d}\rangle_{M}$, the temperature derived from integrated SED (using both single- and multi-temperature model) has higher temperature than $\langle T_{\rm d}\rangle_{M}$ by $\gtrsim 1$ K (except for integrated multi-temperature in the SMC). This difference is because the integrated SED convolves together many regions so it is biased toward warm regions (for single-temperature model), and the limitation in modeling the integrated SED to recover the true distribution of $dM_{d}/dU$ as seen in the highly resolved map (for multi-temperature model).

For the total dust mass, all three methods (integrated SED with single-temperature, integrated SED with multi-temperature, and by summing mass from all pixels in the finest resolution map) agree remarkably well in M31. In the Magellanic Clouds and M33, the mass from the finest resolution map is in between results from the single- and multi-temperature methods. In the Magellanic Clouds, the single-temperature method gives better agreement (but with lower mass) compared to the total mass from resolved map, while in M33, we find better agreement (but with higher mass) between the total masses from multi-temperature method and the resolved map.

7 Summary

Using the archival Herschel maps covering from $\lambda=100$ to $500$ µm and a fitting algorithm developed by Gordon et al. (2014) and implemented by Chiang et al. (2018), we derive maps of equilibrium dust temperature ($T_{d}$), dust mass surface density ($\Sigma_{d}$), and equilibrium infrared luminosity ($\mathcal{J_{\rm IR}}$) for four Local Group galaxies: the Small and Large Magellanic Cloud, M31, and M33. We show these maps, which we make publicly available, in Figure 1. We construct the maps at $13$ pc resolution, the common physical resolution available for the Magellanic Clouds, and $167$ pc resolution, the common physical resolution available for all four targets.

We use these maps to measure how the dust mass and luminosity are distributed as functions of $T_{d}$ (or equivalently, the interstellar radiation field strength, $U\propto T_{d}^{5.8}$), and $\Sigma_{d}$, in each target. We show these distributions in Figures 2 and 3.

Based on these calculations, we gauge the median dust temperatures and surface densities by mass and luminosity and note the characteristic widths of these distributions. We report these values in Table 3. We highlight the following key points.

At $167$ pc resolution, the distribution of dust mass as a function of radiation field implies a median $T_{d}\approx 16{-}20$ K. If we instead consider the distribution of luminosity, the median $T_{d}$ appears higher, $T_{d}\approx 18{-}21$ K. On average, the peak of $T_{d}$ weighting by luminosity is only $\sim 1$ K higher than the peak of $T_{d}$ weighting by mass, but in M33, this shift is about $6$ K, reflecting the strong radial structure in that galaxy. 2. 2.

At $167$ pc resolution, $68\%$ of the dust mass is spread over $\sim 0.4{-}1$ dex in radiation field, $U$. The luminosity spans a lower range, $\sim 0.4{-}0.8$ dex. Again, the strong radial structure and extended disk in M33 lead to a wider distribution, $\sim 0.8{-}1$ dex, than we see in the other galaxies ($\sim 0.4{-}0.6$ dex). 3. 3.

Also at $167$ pc resolution, the dust mass as a function of surface density (i.e., the dust mass surface density PDF) has median value $\log_{10}\Sigma_{\rm d}\approx-1.34$ to $-1.65$, i.e., $\Sigma_{d}\approx 0.02{-}0.05$ M⊙ pc*-2*. This agrees well with observation that the neutral ISM in all of our targets is mostly Hi with dust to gas ratios of a few hundred to a thousand. The median $\Sigma_{d}$ for the distribution of luminosity is usually $\sim 0.1$ dex lower than the distribution by mass. 4. 4.

The widths of the distributions of both mass and luminosity as a function of $\Sigma_{d}$ are $\approx 0.5{-}0.7$ dex. This agrees well with the width measured for gas (atomic and molecular) column density PDFs (e.g., Wada & Norman, 2007; Berkhuijsen & Fletcher, 2015; Corbelli et al., 2018; Sun et al., 2018).

Both the dust mass and luminosity PDFs show significant structure. This structure relates to features visible in the temperature, luminosity, and surface density maps. For example, M33’s strong radial temperature gradient leads to a wide, multi-component distribution in dust mass as a function of temperature. Similarly, M31’s star-forming ring and hot but low-surface density bulge region manifest as clear features in the distribution.

To further quantify how regional features contribute to the global distribution, we identify the active star-forming part of each galaxy and construct a separate distribution for this region and the quiescent region (Figures 4 and 5). We find the following.

As expected, the star-forming parts of galaxies have warmer dust temperature and higher dust mass surface density than the rest of the regions in the galaxy. In our sample, the notable exception is M31, where dust in the central region is heated by older stellar populations (Groves et al., 2012; Smith et al., 2012). This dust stands out in our sample because of its combination of low mass surface density but high dust temperature. 2. 6.

Although these star-forming regions contribute a small fraction to the total dust mass ($10-28\%$), they produce a significant fraction of the equilibrium luminosity in the infrared ($19-51\%$).

Furthermore, we see a more general form of these results by directly plotting the correlation between $\Sigma_{d}$ and $T_{d}$ for the points with high $S/N$ (Figure 6). We find the following.

Dust temperature ($T_{d}$) correlates with dust mass surface density ($\Sigma_{d}$) in SMC, LMC, and M31 (Spearman rank correlation coefficients are weak, but highly significant). In M31, a weak correlation is evident in the dust associated with the star-forming ring, but the hot, low surface density bulge dust heated by the old stellar population leads to an overall anti-correlation between $T_{d}$ and $\Sigma_{d}$. 2. 8.

Very roughly, the sense of the correlation between $\Sigma_{d}$ and $T_{d}$ agrees with expectations from an approximately quadratic gas-star formation scaling relation ($\Sigma_{\rm SFR}\propto\Sigma_{\rm gas}^{2}$), as might be roughly expected in the Hi-dominated regime appropriate for our targets.

The Herschel Space Observatory observed the Local Group galaxies at high physical resolutions. We compare our results to those expected for the same galaxies at farther distances (higher redshift) by progressively convolving the galaxies to coarser resolutions. We find the following.

The luminosity-weighted mean temperature typically exceeds the mass-weighted mean temperature, $\langle T_{\rm d}\rangle_{M}$, by $\sim 2{-}3$ K, an increase of about 10% (or about a factor of $\sim 2$ in $U$). Progressively blurring the maps with a luminosity weighting causes the low resolution $\langle T_{\rm d}\rangle_{M}$ to increase, approaching this higher value (Figure 9). 2. 10.

The best-fit temperature derived from the integrated SED (equivalent to unresolved source) is higher than $\langle T_{\rm d}\rangle_{M}$ from the finest resolution map (Table 4). This bias is somewhat reduced when a sub-resolution multi-temperature distribution (e.g., the approach in Dale et al., 2001; Draine & Li, 2007) is used. 3. 11.

Both points above quantify how severely the temperature from low resolution IR-SED fitting in distant galaxies may overestimate their true mass-weighted mean temperature. As a consequence, the inferred total dust mass in distant galaxies would be underestimated. This is particularly important when the gas mass is determined based on the dust mass.

Lastly, we compare our five-bands fitting results to those that would be obtained by using only one or two bands, as is now commonly required with ALMA (Scoville et al., 2014). We test how the 500 µm intensity, $I_{500\micron}$, performs as a tracer for dust mass surface density, $\Sigma_{d}$. We also test how well our best-fit $T_{d}$ matches individual infrared color ($I_{100\micron}/I_{250\micron}$ or $I_{100\micron}/I_{350\micron}$). We find the following.

A median radio of $I_{500\micron}$ over $\Sigma_{d}$ can recover the trend between $I_{500\micron}$ and $\Sigma_{d}$ in our data (Figure 7). There are two subtleties. First, there is $\lesssim 50\%$ offset in this median ratio from galaxy-to-galaxy. Second, real temperature variations, e.g., as seen in M33 gradient, lead to systematic bias and a modestly sublinear slope between $\Sigma_{d}$ and $T_{d}$. The sense is that without accounting for radial trend in dust temperature, one would estimate that dust in M33 to be more confined to the galaxy center than it actually is. 2. 12.

Although it is an implication of our modeling method to determine the dust temperature, the modified black-body model with a single-temperature and a constant emissivity index, $\beta$, can adequately relate the infrared color and dust temperature (Figure 8), provided that the wavelengths (used as the IR-color) bracket the peak wavelength of the cold dust emission. We observe scatter in this color–temperature relation, particularly in the low temperature. Also, the color–temperature relation in the SMC looks steeper than the prediction from the model, probably due to variation of $\beta$ from region-to-region (Gordon et al., 2014).

We make the maps of dust temperature and mass surface density in the FITS file to be publicly available at https://www.asc.ohio-state.edu/astronomy/dustmaps.

We thank the referee for thoughtful comments that improved the paper. DU, IC, AKL, KMS, and JC are supported by the National Science Foundation (NSF) Grant No. 1615728 and NASA ADAP grants NNX16AF48G and NNX17AF39G. DU and AKL are also partially supported by NSF under Grants No. 1615105, 1615109, and 1653300. Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This work is based, in part, on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Partial support was also provided to AKL in the context of HST project HST-GO-13659.013-A (PI Sandstrom), which is studying dust in the Magellanic Clouds. That project is based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555.

Appendix A Challenge in comparing across physical resolutions

One challenge in fitting the dust properties across physical resolutions is computing the covariance matrix of the background fluxes ($\S$3.2) for the low-resolution images. The convolution process blurs the image. Therefore, to conserve flux from a galaxy, the background area decreases and the region of interest area increases in the low resolution maps. This means we can only compute the covariance matrix of background flux correctly down to certain resolution. In this study, this lower limit in resolution is selected to be where there are at least ten resolution elements across the target image.

It is shown in Chiang et al. (2018) that the choice of covariance matrix heavily affects the results of fitting. Thus without the covariance matrix of background fluxes, we are unable to fit dust properties at coarser resolutions. Several attempts have been made to extrapolate the covariance matrix to coarser resolutions. So far, none of the methods we try are proved to give the correct matrix mainly due to two reasons. (1) The background flux consists of both random noise and real signal from background objects. (2) The observed correlations between Herschel bands are not constant across resolutions. Whether it is possible to extrapolate the covariance matrix for low-resolution images from high-resolution images remains an open question for future study.