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

, where β is the dust emissivity index. The temperature-ISRF relation assumes that the dust cross section per unit mass, κ, depends on the wavelength as κ∝λ -β (Draine 2011).

Extragalactic observations necessarily convolve together many different environments and radiation fields because of the limited angular resolution of infrared (IR) telescopes. To account for this, starting from the work of Dale et al. (2001), many extragalactic studies have modeled the IR 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 d /dU, and is usually modeled as a power law, U -α , where α∼2 (Draine et al. 2007).

Some part of this distribution of dust heating, dM d /dU, reflects geometry and radiative transfer on very small scales. For example, in the subparsec 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.

Because of 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 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 objects, which mimics the high-redshift studies (e.g., Magdis et al. 2010Magdis et al. , 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), and Chastenet et al. (2017); the M31 work by Fritz et al. (2012), Groves et al. (2012), Smith et al. (2012), and Draine et al. (2014); and the M33-focused investigations of Braine et al. (2010), Boquien et al. (2011), and 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 (∼10-100 pc) structures for observations of more distant galaxies, and employing the new fitting code developed by Chiang et al. (2018), which is based on Gordon et al. (2014).

Throughout this paper, we use the five longest wavelengths of Herschel bands (100, 160, 250, 350, and 500 μm). We expect that the IR 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 Section 2 and our modified blackbody modeling in Section 3. We present maps of dust mass surface density and temperature for all four targets in Section 4. We also show the distribution of dust mass as functions of illuminating radiation field for each galaxy and locations within a galaxy in Section 4. We show the correlation between dust mass surface density versus dust temperature, IR color versus dust temperature, and 500 μm intensity versus dust mass surface density in Section 5. We explore how the inferred dust temperature changes as a function of resolution (down to treating the galaxy as one unresolved source) in Section 6. Finally, we summarize our findings in Section 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°. 6 for SMC (Subramanian & Subramaniam 2012), 34°.7 for LMC (van der Marel & Cioni 2001), 77°. 7 for M31 (Corbelli et al. 2010), and 56°for M33 (Paturel et al. 2003). We summarize the symbols used in this paper in Table 1.

Region of Interest

We identify a region of interest for each galaxy on the basis of 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 fit only 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 three times the noise level as inferred from the 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 ratio (S/N) cut of 2 or 4 and also by 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 varies within only few percent by using a smaller or bigger mask. The region of interest definition has a greater 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.

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 IR background and resolved background galaxies, and the adopted observing and imaging strategies that 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 H I 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 drop pixels that deviate by more than ±3σ from the fit, where σ is the robustly estimated rms noise. Exclusion of these outlier pixels is important to exclude pixels that may be contaminated by instrumental artifacts or nonbackground 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 represents a small adjustment (significantly less than 1σ).

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

Convolution and Reprojection

We convolve the background-subtracted 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 2000 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 , and dust mass surface density (right panels) at 13 pc resolution for both LMC and SMC, and at 167 pc resolution for both M31 and M33. White contours mark the region of interest (Section 2.2). In the LMC, SMC, and M33, high dust temperatures appear associated with star-forming regions. In M31, the old stellar populations heats the dust in the bulge to higher temperatures (Groves et al. 2012;Draine et al. 2014;Viaene et al. 2014). In all targets, the dust surface density maps visually resemble other maps of the interstellar matter (ISM), especially atomic hydrogen (see Kim et al. 1998;Stanimirovic et al. 1999;Staveley-Smith et al. 2003;Braun et al. 2009;Braun 2012;Koch et al. 2018).

beam (i.e., we oversample at roughly the Nyquist sampling rate). This means the pixel size in parsecs 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).

Modified Blackbody Model

The intensity of dust emission, I ν , for dust in equilibrium with the radiation field depends on the dust optical depth, τ(λ), and the dust equilibrium temperature, T d , via

where B ν (T d , λ) is the Planck function for temperature T d at a wavelength λ.

For our target galaxies and resolutions, dust emission at λ100 μm is almost always optically thin. In this case,

The dust optical depth can be written as

where κ(λ) is the dust absorption cross section per unit mass at a wavelength λ, and Σ d is the dust mass surface density.

Fitting the SED yields τ(λ), but for a given τ(λ), the values of κ(λ) and Σ d are degenerate. Breaking this degeneracy requires an independent calibration of κ(λ). Such a calibration in turn requires that τ(λ) be measured in a location where Σ 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), which is based on Gordon et al. (2014). We calculate τ(λ) by fitting the solar neighborhood cirrus emission using our modified blackbody model. Then, we derive κ(λ) from τ(λ) by adopting the dust abundance found by Jenkins (2009). We further assume that κ(λ) remains constant across all regions in Local Group galaxies. Though variations in κ do represent a systematic uncertainty, this procedure removes another systematic uncertainty by using the same dust model to estimate κ 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

where λ 0 is a reference wavelength, and β represents the dust emissivity index. In this paper, we assume a fixed β=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 Σ d from different choices of β values would be small because κ(λ 0 ) is calibrated accordingly.

Combining those assumptions above, we get

We adopt λ 0 =160 μm as our reference wavelength. For β=1.8, our fit to the Milky Way cirrus yields κ(λ 0 )= 18.7±0.6 cm 2 g -1 .

Fitting the Dust SED

For each line of sight, we estimate Σ d and T d by simultaneously fitting the intensities at 100, 160, 250, 350, 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, , given the observed SED. Here,

where ΔI indicates a vector recording the difference between the observed SED and the model SED. Before calculating ΔI, we integrate the models (Equation ( 5)) with the appropriate response function of that Herschel band.  -1 refers to the the inverse of the covariance matrix. This term takes into account the interband 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, , is calculated by summing the covariance matrix calculated 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, , for each model across a grid of T d and Σ d . The grid for T d spans from 5 to 50 K with increments of 0.5 K. We space the grid for Σ d logarithmically, covering from -4 to 1 dex in log 10 Σ d with a step size of 0.025 dex.

From the grid of likelihoods, we calculate the expectation values for T d and Σ d via

Here, the sum covers all cells in the fitting grid. We adopt these expectation values as our estimate of the best-fit Σ d and T d at each pixel inside the region of interest. We correct Σ d for galaxy inclination, i, by multiplying Σ 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.

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 κ, T d relates to the ISRF, U,

The normalization takes T d =18 K to corresponds to the solar neighborhood radiation field, U=1 (Draine et al. 2014). This normalization agrees well with our fit to calibrate κ. In that case, we found T d =18.29±0.11 K for the local cirrus.

Assuming that the solar neighborhood radiation field heats local cirrus, then this supports equating T d ≈18 K with U≈1.

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

Here,  IR represents the IR luminosity surface density of grains with mass surface density Σ d at temperature T d . We refer to  IR as the "equilibrium luminosity surface density." Even though the IR luminosity captures the integral of IR SED under our best-fit model, it underestimates the true total IR 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 make only a small difference in the total dust mass. However, they produce an important fraction of the total IR luminosity emitted by all kind of dust.

Caveats

We adopt several simplifying assumptions that the reader should bear in mind. We assume a single T d for each pixel. In reality, T d may vary within any given resolution element, especially in star-forming regions, which will contain many local heating sources. Hence, our T d measured at 13 or 167 pc resolution is a luminosity-weighted mean of the subresolution T 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 d (Draine & Li 2007) are beyond the scope of our work and will be addressed in upcoming paper (J. Chastenet et al. 2019, 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 β also affect the best-fit T d , with the sense that a higher value of β gives a lower value of the best-fit T d (Kelly et al. 2012). There are some measurements suggesting a variation of β within and among Local Group galaxies (e.g., Smith et al. 2012;Gordon et al. 2014;Tibbs et al. 2018). We fix β in the interests of simplicity and to avoid the ambiguities in 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 κ(λ 0 )=18.7±0.6 cm 2 g -1 for λ 0 =160 μm and β=1.8. This value relies on the assumption that the optical properties of dust grains throughout the Local Group galaxies are the same as those in the Milky Way cirrus. If κ(λ 0 ) is underestimated, we will overestimate the dust mass, and viceversa. 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 κ(λ 0 ) varies from galaxy to galaxy, but not between regions within a galaxy, our results will still be robust.

Maps and Distributions of Dust Mass, Luminosity, and Temperature

Figure 1 shows our best estimates of dust equilibrium temperature, equilibrium IR luminosity surface density, and dust mass surface density for each target. 3 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 (see 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 starforming 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., Gordon et al. 2006;Nieten et al. 2006), and in M33, the dust mass maps show flocculent spiral structures.

Atomic hydrogen, H I, makes up most of the neutral ISM in each of our targets. As a result, our Σ d maps resemble H I 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 to date.

The luminosity maps ( º S T d d IR 5.8 ) show the brightest regions where most IR light comes from. These bright regions are due to high Σ 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 d and dense Σ d . M31 shows distinctive feature, where the bulge (that shows a "hole" in the Σ d map) is luminous because of its high T d .

On the basis of these maps, we calculate the distributions of dust mass and luminosity as functions of local illuminating radiation field, µ U T d 5.8 (in Section 4.1), and local dust mass surface density, Σ d (in Section 4.2). We compare those distributions in star-forming regions against the rest of area in galaxies in Section 4.3.

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

The dust mass distribution as a function of illuminating ISRF strongly affects the integrated SED of a galaxy. Though we cannot measure the subresolution scale distribution (e.g., photodissociation region or H II region structure), our maps allow us to measure the dust mass distribution as a function of illuminating radiation field, 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 d /dlog 10 U, where dM d is the amount of dust mass heated by U within a range of ±0.5dlog 10 U. We choose dlog 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 d by the total dust mass, M Tot , in each galaxy, and define the mass fraction as dm d =dM d /M Tot .

Then, the normalized distribution function is

We interpret this normalized distribution function, dm d /dlog 10 U, as the fraction of the total dust mass per dex of U, calculated in a bin with width of dlog 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 and3, which shows the 16th-to-84th percentile across all realizations. The Figure 2. Distributions of dust mass (left) and equilibrium IR luminosity (right) at 13 pc resolution for the Small and Large Magellanic Clouds as functions of (top) illuminating radiation field, U, and (bottom) dust mass surface density, Σ d . We calculate the distributions in 0.01 dex-wide bins, normalized to the total mass or total luminosity in the galaxy, and then smooth the distribution using a five bin-wide Gaussian filter. Shaded regions indicate 16th-to-84th percentile confidence. Circles, ticks, and lines at the top of the plot the peak, median, and 16th-to-84th percentile width of the distribution. Characteristic widths for all four distributions are ∼0.6-0.95 dex, but the detailed distributions differ. The luminosity distributions tend to reflect higher-mass surface density and higher radiation field regions than the mass distribution. Both galaxies show significant contributions from individual bright regions to the luminosity distribution; the "bumps" in the distribution reflect the bight spots in the luminosity maps (middle column) of Figure 1. We quantify these distributions in Table 3. We compare the distributions in the Magellanic Clouds to those in M31 and M33 at the common resolution of 167 pc in Figure 3.

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 Σ 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 Σ 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 ∼±1σ uncertainty as the difference between the 84th and 16th percentile value across all realizations in that bin.

Parameters of the Distribution. From each distribution, we measure three quantities described as follows.

1. U peak , defined as the value of U that corresponds to the peak of the distributions of mass or luminosity. This is the most common value, by mass or luminosity, in the distribution function. It often, but not always, lies near the center of the distribution. 2. U med , defined as the median value of U in the distributions of mass or luminosity. 3.  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 peak and U med for each galaxy by a circle and a tick above the distributions in Figures 2 and3, and the width as horizontal lines, color-coded by galaxy. We record their values in Table 3.

Distribution of Dust Mass as a Function of Radiation Field

In practice, we calculate dm d /dlog 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 Figure 3. Distributions of dust mass (left) and equilibrium IR luminosity (right) at 167 pc resolution for all four targets. The plots show distributions of dust mass (left) and equilibrium IR luminosity (right) as functions of (top) illuminating radiation field, U, and (bottom) dust mass surface density, Σ d . We calculate the distributions in 0.01 dex-wide bins, normalized to the total mass or luminosity in the galaxy, and then smooth the distribution using a five bin-wide Gaussian filter. Shaded regions indicate 16th-to-84th percentile confidence. For the SMC and LMC, these are larger than in Figure 2 because of the smaller number of pixels at this resolution. Circles, ticks, and lines at the top of the plot show the peak, median, and 16th-to-84th percentile width of the distribution, respectively. The majority of dust mass and luminosity in Magellanic clouds originates from regions with ISRF higher than the solar neighborhood value (log 10 U e =0). In M31 and M33, most of the dust mass has U<U e , but this is different in the luminosity distribution, where most luminosity in M31 at around U e and higher than U e in M33. mass in it, we assign the value of dm d /dlog 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 and3.

Because our estimate of U depends directly on our best fit of dust temperature, T d , via Equation (10), our measured distribution function can also be converted to the distribution of mass as a function of

In the top left panels of Figures 2 and3, 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 peak ≈-0.01 dex (corresponds to T peak ≈18 K). This is very close to the solar neighborhood value of log 10 U e =0 dex. However, as the top left plot in Figure 3 shows, U peak in M31 lies at the upper end of a wide distribution (i.e., U peak >U med ≈-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 peak in our samplefoot_2 (-0.52 dex), but this peak lies toward the low end of a wide distribution (i.e., U peak <U med ≈-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 toward the inner galaxy than the dust mass.

At 167 pc resolution, both Magellanic Clouds show narrower distributions and higher U peak than M31 and M33. The LMC has the highest log 10 U peak (0.30 dex) in our sample, while the SMC has log 10 U peak of 0.16 dex. These values are close to their log 10 U 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=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 Section 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).

Distribution of Luminosity as a Function of Radiation Field

We also calculate the distribution of equilibrium IR luminosity as a function of U, d IR /d log 10 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,  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 and3.

As one might expect, given the dependence of  IR on µ U T d 5.8 (Equation ( 11)), the luminosity distributions shift toward higher U compared to the mass distributions. This strongly affects U peak in M33, which has the lowest U peak and U med by mass but has the highest U peak by luminosity (0.44 dex). The other galaxies show modest shifts in U peak and U med , but the shape of their luminosity distributions around the peak change significantly (compared to that in the mass distributions). 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 Section 4.3).

At 167 pc resolution, M31 has the lowest log 10 U peak (0.04 dex) and lowest log 10 U 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 med for luminosity is -0.03, but only -0.14 for mass distribution (i.e., lower value of U 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, in Figure 1, the central part of the distribution, which maps to the star-forming 10 kpc ring (described in Section 4.3), 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%).

M33 shows the most striking change between the mass and luminosity distributions. Its log 10 U peak shifts from the lowest (-0.52) by mass to the highest (0.44) by luminosity. Its U 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 d , along with the extended gas and dust surface density distributions (e.g., Koch et al. 2018). As the maps (Figure 1) show, the mass in M33 extends out to large radii, while the luminosity appears more centrally concentrated, coincident with active star formation (Section 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 starforming regions (high-end tail in the top right panel of Figure 2). Even though warm dust components can contribute toward these hot regions, they contribute only a negligible fraction in the mass distribution (top left panel of Figure 2) and lies outside 68% of the mass distribution. Hence, the accuracy of our modeling for these warm dust component does not affect the main results of this paper.

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-todust ratio (GDR). 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 Σ d ∼0.01-0.1 M e pc -2 . For a Galactic GDR of ∼150, this range equates to a gas mass surface density of ∼1-15 M e pc -2 . The metallicities of the LMC, SMC, and M33 are all lower than that of 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 ∼3-50 M e 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., Wong et al. 2011;Druard et al. 2014).

The LMC and M33 show similar mass distributions. They exhibit the highest Σ peak of ≈0.045 M e pc -2 in our sample (their Σ peak ≈Σ med ) and widths of 0.5-0.7 dex. These two dwarf spirals have comparable metallicity, Z≈0.5Z e (Russell & Dopita 1992;Rosolowsky & Simon 2008), and the distribution appears consistent with dust being mostly mixed with high column density of H I in both targets.

The SMC shows low Σ peak and low Σ med of 0.02 and 0.03

M pc 2 , respectively. The SMC has a lower metallicity and higher GDR of ∼1200, compared to those of M33 and the LMC (≈380; Roman-Duval et al. 2014). The lower Σ peak in the SMC could thus be expected if the gas mass surface density in SMC were comparable to that in the LMC. In reality, the SMC's elongated structure along the line of sight leads to higher H I column densities than one finds in the other Local Group targets (Stanimirovic et al. 1999;Staveley-Smith et al. 2003). These effects combine to produce the multicomponent distribution 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 Σ peak and Σ med ∼0.02 M e 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°.7; Corbelli et al. 2010).

Implications

These mass and luminosity distributions as a function of Σ d have two implications: low IR optical depth and no evidence for opaque H I, as described below.

Low IR Optical Depth. The dust mass surface densities in our targets almost never exceed Σ d ∼1 M e pc -2 . Given our adopted value of κ at 160 μm, this implies that at our resolutions, the optical depth throughout our targets at 100 and 160 μm will be 1×10 -2 and 5×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 IR opacity, τ IR (e.g., see Thompson et al. 2005), and should exceed the UV-driven radiation pressure when τ IR >1. Thus, we expect radiation pressure from long-wavelength IR photons will contribute an insignificant amount of feedback over essentially all of the area in these targets.

The highest Σ d that we observe in the LMC (Figure 2) implies significant opacity at shorter IR wavelengths (e.g., λ∼20 μm). We do not consider emission at these wavelengths. 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 IR photons in the hottest, highest-density parts of our targets (Lopez et al. 2011(Lopez et al. , 2014)). To see this, note that at 13 pc resolution, the LMC by luminosity includes a small, but noticeable, contribution at Σ d ∼1-10 M e pc -2 (bottom right panel of Figure 2). This Σ d ∼1 M e pc -2 implies τ 160 μm of only ∼4×10 -3 . If we extrapolate to ∼20 μm, though, τ 20 μm ∼0.15-0.3 (the extrapolation using β∼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 nontrivial, and perhaps even dominant, contribution to the radiation pressure from IR-reprocessed photons.

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

In contrast, our Σ d maps do not have a similar patchy appearance, nor do the values of Σ 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 GDRs of ∼100-500 (Leroy et al. 2011;Roman-Duval et al. 2014). In those cases, the Σ d that we find corresponds to column densities of <5×10 22 cm -2 . In order for the opaque H I features posited by Braun (2012) to exist but not appear in our Σ d maps, they would need to be unusually dust poor (i.e., GDR >1000). 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 GDR or the presence of COdark molecular gas. This issue certainly remains an open topic. Our maps simply do not support with specific structure found by Braunʼs analysis.

Regional Comparisons

As we see in Figures 2 and3, 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 Sections 4.1 and 4.2). For M31, we also build separate distributions for the hot, low Σ d central region (inside 1 kpc from the nucleus).

Selection of Regions. We define "star-forming" regions as areas with Spitzer MIPS (Rieke et al. 2004;Werner et al. 2004) 24 μm intensity >1 MJy sr -1 (shown as the blue regions in the right panels of Figures 4 and5). This intensity is equivalent to » ´- M 7 10 4 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. 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°.7 from Corbelli et al. 2010).

Resulting Distributions. We show the mass and luminosity distribution functions of star-forming regions (blue) and the rest of area in the galaxy (red) in Figures 4 and5 (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 multicomponent distributions in Figures 2 and3 are a result of superimposing physically distinct regions. Star-forming regions, which also tend to be the regions rich in molecular gas, contribute many of the high-Σ d and high-U features in the full-galaxy distributions. Because the star-forming regions have both high T d and high Σ d , they appear to be 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 starforming ring or the rest of the galaxy (≈30 K versus 17 K), but it captures only a very small fraction of the total dust mass (≈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.

Correlations between S d and T d

Both the maps (Figure 1) and the distributions (Figures 2345) suggest a correlation between Σ 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 S/N, we expect correlated uncertainties to drive an apparent anticorrelation between the best-fit Σ 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 S/N less than 3 in either T d or Σ 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×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 Σ d , in regions where log 10 Σ 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 that warm regions also tend to have high mass surface densities. As we saw in the last section, starforming regions in the maps tend to have both high mass surface density and high temperature. To some degree, a correlation between T d and Σ d should be expected on large scales for dust internally heated by star formation. Following Schmidt (1959) and Kennicutt (1998) 5.8 and Σ gas ∝Σ d , then the correlation between Σ d and T d 5.8 appears when the Σ SFR -Σ gas relationship is steeper than linear. Indeed, these galaxies are all mostly H I-dominated, and the Σ SFR -Σ H I scaling relation does have a superlinear 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 Σ d and T d 5.8 is shown as the black line in Figure 6 (with an arbitrary normalization).

This S µ T d 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 because of 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 M pc 2 for dust surface density. The gray contours mark the probability of data points of 0.01, 0.1, and 0.5 points per grid, while the blue contours mark the probability density data points of 0.01, 0.1, and 1.0 points per grid. The red points show the median temperature in 0.1 dex bin in log 10 Σ d >-1.5. The solid black line show a relation of 5.8log 10 T d ≈9+log 10 Σ d . The Spearman (1904) rank correlation coefficient, r, for all blue points is indicated. In the SMC, LMC, and M33, the temperature and mass for blue contours are correlated, where high mass surface density tends to be located in the warm star-forming regions. However, in M31, dust temperature is uncorrelated with the dust mass surface density.

discrete features. This effect can also be seen in the Magellanic Clouds (Jameson et al. 2016). In this case, the local Σ d will no longer be directly correlated to the amount of heating.

Intensity at 500 μm and Dust Mass Surface Density

With powerful submillimeter telescopes like Atacama Large Millimeter Array (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 Σ d from our five bands fitting, which leverages information from 100 to 500 μm SED, to the monochromatic intensity at 500 μm. This offers a check on the reliability of submillimeter intensities as a dust mass tracer.

Figure 7 shows the correlation between resolved I 500 μm and Σ d over three orders of magnitude in both axes. We select only pixels above 2σ uncertainty for both I 500 μm and Σ d . Contours show the density of data points, generated using a Gaussian kernel density estimate with a kernel size of ≈0.02 dex. In each panel, lines show relations with a power-law index of 1 with normalization set by the median of Σ d /I 500 μm measured for each galaxy. The value of this median ratio is Σ d /I 500 μm ≈ - I 500 μm 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 to linearly trace Σ 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 predicts Σ 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 GDR.

IR Color and Temperature

We compare the 100-to-250 μm or 100-to-350 μm IR 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 μm /I 250 μm and I 100 μm /I 350 μm as IR colors because those wavelengths bracket the peak of the dust SED at ∼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 Section 5.2. We also plot the relation between temperature and colors expected from modified blackbody models with β=1.4, 1.8, and 2.2.

The data broadly follow the model, though with some notable deviations. The plots suggest that LMC and M33 may be better described by β closer to ∼2 than our adopted β∼1.8. Both the SMC and the LMC show points covering a range of β, with the regions with redder color more consistent with higher β, while the rest of the pixels are in good agreement with our adopted value of β. Note that despite the shift in β, these features are as likely to 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 β). 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).

We note that in the range of 30<T d <40 K, the modified blackbody models still show a curvature in 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 (dotted), 1.8 (solid), and 2.2 (dashed). The Magellanic clouds are at 13 pc resolution, while M31 and M33 are at 167 pc resolution. M31 data follow β=1.8 model very well but others show deviations. LMC and M33 look better to be described with 1.8<β<2.2, while SMC shows variation in β with warm regions have β∼1.8 and cold regions have β∼2.2.

tail), and this is clearly not the case in our study. Therefore, the high-temperature tail in the distributions (Figure 2) is real.

Dust Temperature as a Function of Physical Scale

We study four of the neighboring galaxies at the diffraction limit of the best IR telescope to date. Such high physical resolution is not available for distant galaxies. Still, as discussed in Section 1, there remains considerable interest in using dust as a tracer of the ISM 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 from convolving the SED first and then fitting that convolved SED. The reasoning behind this choice is that we lose the background area gradually in the coarser resolution maps. This would make the covariance matrix calculation (Equation ( 7)) difficult, which leads to nonrobustness of the fitting results of the convolved SED (see Appendix).

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-weighted (rather than mass-weighted) averages at some coarser resolution. 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, Σ d,L , at resolutions coarser than our original maps. To calculate T d,L , we begin with our highest resolution T 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

where  IR is a quantity proportional to the equilibrium dust luminosity per unit area in the highest resolution map (as defined by Equation (11)), * denotes convolution, and  is the Gaussian kernel to convolve from the finest resolution to the target coarser resolution.

Similarly, we calculate the luminosity-weighted mass surface density, Σ d,L , via

In Equations ( 14) and (15), we calculate T d,L and Σ 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

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

Here, the sum å i runs over all pixels inside the region of interest.

As we see in Section 4.3, weighting by luminosity leads to higher á ñ T d 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 "subresolution" 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 á ñ T d M (black dots) for each galaxy as a function of physical resolution. We also mark á ñ

T d M at the finest resolution as blue line. 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. Because of the different physical sizes among galaxies, LMC and SMC are convolved to seven resolutions from 13 to 500 pc, while M31 and M33 are convolved to five resolutions from 167 pc to 2 kpc.

In general, á ñ T d 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 á ñ T d M toward higher temperatures. We already saw an indication of this when comparing the distributions of mass by radiation field at 13 and 167 pc resolution for the Magellanic Clouds (Figures 2 and3).

This trend with spatial resolution is weakest in M31; the inferred mean temperature changes by only <1 K as we blur the maps by a factor of ∼10 in resolution. We identify two likely causes. First, we see a weaker correlation between temperature and mass surface density in Figure 6. Second, much of the variation in temperature occurs radially on scales larger than a kiloparsec, so that blurring together the distinct physical conditions is less of an issue in M31.

The SMC and M33 show the largest increases in á ñ T d M with coarser spatial resolution. The implied temperature in both galaxies increases by ∼2 K as the resolution degrades by a factor of ∼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, á ñ T d L (red), and massweighted temperature, á ñ T d M (blue), at the finest resolution. At the finest resolution, the black point lies on the blue line by construction. As we blur the galaxy toward 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 ∼2-3 K, or about 10%, change in T d , or a factor of ∼2 change in U.

Implication for the Total Mass Measurement. Luminosity weighting tends to bias toward 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 ν ∝T d Σ d . As the luminosity weighting biases toward higher temperature, it will also lead to an underestimate of the mass.

To a first order, the mass will be biased low by the same fraction as the temperature is biased high. This is a ∼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.

For the LMC, Galliano et al. (2011) considered this effect and estimated that the mass inferred from the integrated SED is biased low by 50%. They included shorter wavelength bands than 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 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 Figure 9. Mass-weighted mean of dust temperature (black dots) as a function of resolutions for different galaxies. We mark the mass-weighted and the luminosityweighted mean temperature measured at the finest resolution as blue and red lines, respectively. The green and orange lines are temperature from fitting the integrated SED with single-and multitemperature approaches, respectively. The error bars and shaded areas are derived by taking the 16th-to-84th percentile range in the probability distribution function of the likelihood. This figure shows that the mass-weighted mean of temperature increases in coarser resolutions and approaching the luminosity-weighted average (red line) because convolution favors bright regions that are also warm.

accurate because a significant portion of blending populations (unresolved by our study) must still occur even below that resolution.

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. 2010Magdis et al. , 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; singletemperature and multitemperature modified blackbodies. Figure 10 shows these integrated SEDs as red points, along with the best-fit models (curves). Uncertainty in the best-fit models are visualized by drawing 100 realizations in proportion to their likelihood and 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 multitemperature blackbody, and (3) the mass-weighted mean of temperature, á ñ T d M , calculated from the finest resolution map.

Single-temperature Model. First, we use the same singletemperature modified blackbody model 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 from the luminosity-weighted mean of temperature, á ñ T d L , because in the integrated SED, we take average of intensity without any weighting. Therefore, we do not expect á ñ T d L to be the same as the temperature derived from the integrated SED fitting (red versus green lines in Figure 9).

The expected values of T 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 that the unresolved observations would overestimate the temperature of the majority of the dust by 3 K or a factor of 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 17% for an object with 18 K temperature.

Multitemperature Model. We also fit the integrated SED by adopting a subresolution multitemperature distribution in the form of (for a ¹ 1; Draine & Li 2007)

where δ is Dirac's delta function. Here, a fraction γ of the dust mass within a resolution element is heated by a power-law distribution field U -α between U min and U max , while the rest is heated by the ambient ISRF of U min . We use the same grid-based fitting procedure as before (Section 3.2). The grid for α is from 1.1 to 5.1 with an increment of 0.2, the grid for log 10 γ is from -3 to 0 with an increment of 0.2, and the grid for log 10 U min is 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 (Section 3.2). The value of U max is fixed at 10 3 . The adoption of this U max value does not affect the outcome of fitting parameters.

The best-fit model is shown as the dashed purple curve in Figure 10. The expected values for U min , α, and γ are listed in Table 5. The value of U min represents the ambient ISRFs that illuminate the majority (1-γ) 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 (Section 4.1). The rest of dust mass (a fraction of γ) is heated by ISRF in the range between U min and U max .

The value of α 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 min <U<U max . A shallower slope (smaller α), 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 min and U max . Hence, this multitemperature 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

Then, we convert á ñ U M to the dust temperature by inverting Equation (10).

Comparison between methods. In SMC and M33, the temperature derived from the integrated SED multitemperature model is lower than those derived from the integrated SED single-temperature model by ∼1-2 K, while the opposite happens in M31. The best-fit temperature from the two methods above matches in LMC. Compared to á ñ T d M , the temperature derived from integrated SED (using both singleand multitemperature model) has higher temperature by 1 K (except for the integrated multitemperature model in the SMC). This difference is because the integrated SED convolves together many regions so it is biased toward warm regions (for single-temperature), 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).

For the total dust mass, all three methods (integrated SED with single-temperature, integrated SED with multi-T d , 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 between the results from the single-and multitemperature methods. In the Magellanic Clouds, the single-temperature method gives better agreement (but with lower mass) compared to the mass from the resolved map, while in M33, we find better agreement (but with higher mass) between the total dust masses from multitemperature method and the resolved map.

61.7±23.3

Note. For multi-T, we use parameters U min , U max , and γ to derive the mean ISRF, á ñ U , using Equation (19), and then convert it to dust temperature using Equation (10).

and implemented by Chiang et al. (2018), we derive maps of equilibrium dust temperature (T d ), dust mass surface density (Σ d ), and equilibrium IR luminosity ( IR ) for four Local Group galaxies: the SMC, LMC, M31, and M33. We show these maps, which we make publicly available, in Figure 1. We construct the maps at 13 pc resolution, which is the common physical resolution available for the Magellanic Clouds, and at 167 pc resolution, which is the common physical resolution available for all four targets. We use these maps to measure how dust mass and luminosity are distributed as functions of ISRF intensity, U (or equivalently T d ), and Σ d , in each target. We show these distributions in Figures 2 and3.

On the basis of 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.

1. At 167 pc resolution, the distribution of dust mass as a function of radiation field implies a median T d ≈16-20 K.

If we instead consider the distribution of luminosity, the median T d appears higher, T d ≈18-21 K. On average, the peak of T d weighting by luminosity is ∼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. At 167 pc resolution, 68% of the dust mass is spread over ∼0.4-1 dex in radiation field, U. The luminosity spans a similar range, ∼0.4-1 dex. Again, the strong radial structure and extended disk in M33 lead to a wider distribution, ∼0.8-1 dex, than we see in the other galaxies (∼0.4-0.6 dex). 3. At 167 pc resolution, the dust mass as a function of dust mass surface density (i.e., the dust mass surface density PDF) has median value of Σ d ≈0.02-0.05 M e pc -2 . This agrees well with observation that the neutral ISM in all of our targets is mostly H I with dust-to-gas ratios of a few hundred to a thousand. The median Σ d for the distribution of luminosity is usually ∼0.1 dex higher than the distribution by mass. 4. The widths of the distributions of both mass and luminosity as a function of Σ d are ≈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 mass surface density maps. For example, M33ʼs strong radial temperature gradient leads to a wide, multicomponent distribution in dust mass as a function of temperature. Similarly, M31ʼs star-forming ring and hot but low-mass 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 and5). We find:

1. As expected, the star-forming parts of galaxies have warmer dust temperature and higher mass surface density than that in the rest of area 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. 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 IR (19%-51%).

Furthermore, we see a more general form of these results by directly plotting the correlation between Σ d and T d for the points with high S/N (Figure 6). We find:

1. Dust temperature (T d ) and mass surface density (Σ d )

positively correlate with a low, but highly significant, Spearman rank correlation coefficient in all targets except M31. 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 anticorrelation. 2. Very roughly, the sense of the correlation between Σ d and T d agrees with expectations from an approximately quadratic gas-star formation scaling relation (Σ SFR ∝Σ gas 2 ), as might be roughly expected in the H I-dominated regime appropriate for our targets.

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 convolving the galaxies to progressively coarser resolutions. We find:

1. The luminosity-weighted mean temperature typically exceeds the mass-weighted mean temperature, á ñ T d M , by ∼2-3 K, an increase of about 10% (or about a factor of ∼2 in U). Progressively blurring the maps with a luminosity weighting causes the low resolution á ñ T d M to increase, approaching this higher value (Figure 9).

2. The best-fit temperature derived from the integrated SED (equivalent to unresolved source) is higher than á ñ T d M from the finest resolution map (Table 4). This bias is somewhat reduced when a subresolution multitemperature distribution (e.g.,the approach in Dale et al. 2001;Draine & Li 2007) is used.

Both results 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 total dust mass would be underestimated. This is particularly important in high-redshift studies where the gas mass is derived based on the dust mass.

Finally, we compare our five band 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 μm , performs as a tracer for dust mass surface density, Σ d . We also test how well our best-fit T d matches individual IR color (I 100 μm /I 250 μm or I 100 μm /I 350 μm ). We find:

1. A median radio of I 500 μm over Σ d can recover the trend between I 500 μm and Σ d in our data (Figure 7). There are two subtleties. First, there is ∼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 Σ d and T d . The sense is that without accounting for dust temperature, one would estimate the dust in M33 to be more confined to the galaxy center than it actually is. 2. Although it is an implication of our modeling method to determine the dust temperature, modified blackbody model a single temperature and a constant β can adequately relate the IR 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 β from region-to-region (Gordon et al. 2014).

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