Modeling Dust and Starlight in Galaxies Observed by Spitzer and Herschel: The KINGFISH Sample

Interstellar dust and starlight are modeled for the galaxies of the project “Key Insights on Nearby Galaxies: A Far-Infrared Survey with Herschel.” The galaxies were observed by the Infrared Array Camera and the Multiband Imaging Photometer for Spitzer on Spitzer Space Telescope, and the Photodetector Array Camera and Spectrometer and the Spectral and Photometric Imaging Receiver on Herschel Space Observatory. With data from 3.6 to 500 μm, dust models are strongly constrained. Using a physical dust model, for each pixel in each galaxy we estimate (1) dust surface density, (2) dust mass fraction in polycyclic aromatic hydrocarbons (PAHs), (3) distribution of starlight intensities heating the dust, (4) total infrared (IR) luminosity emitted by the dust, and (5) IR luminosity originating in subregions with high starlight intensity. The dust models successfully reproduce the observed global and resolved spectral energy distributions. With the angular resolution of Herschel, we obtain well-resolved maps (available online) for the dust properties. As in previous studies, we find the PAH fraction to be an increasing function of metallicity, with a threshold oxygen abundance Z/Z⊙ ≈ 0.1, but we find the data to be fitted best with increasing linearly with above a threshold value of 0.15(O/H)⊙. We obtain total dust masses for each galaxy by summing the dust mass over the individual map pixels; these “resolved” dust masses are consistent with the masses inferred from a model fit to the global photometry. The global dust-to-gas ratios obtained from this study are found to correlate with galaxy metallicities. Systems with Z/Z⊙ ≳ 0.5 have most of their refractory elements locked up in dust, whereas in systems with Z/Z⊙ ≲ 0.3 most of these elements tend to remain in the gas phase. Within galaxies, we find that is suppressed in regions with unusually warm dust with . With knowledge of one long-wavelength flux density ratio (e.g., f160/f500), the minimum starlight intensity heating the dust ( ) can be estimated to within ∼50%, despite a variation in of more than two orders of magnitude. For the adopted dust model, dust masses can be estimated to within ∼0.2 dex accuracy using the f160/f500 flux ratio and the integrated dust luminosity, and to ∼0.07 dex accuracy using the 500 μm luminosity alone. There are additional systematic errors arising from the choice of dust model, but these are hard to estimate. These calibrated prescriptions for estimating starlight heating intensity and dust mass may be useful for studies of high-redshift galaxies.


Abstract
Interstellar dust and starlight are modeled for the galaxies of the project "Key Insights on Nearby Galaxies: A Far-Infrared Survey with Herschel." The galaxies were observed by the Infrared Array Camera and the Multiband Imaging Photometer for Spitzer on Spitzer Space Telescope, and the Photodetector Array Camera and Spectrometer and the Spectral and Photometric Imaging Receiver on Herschel Space Observatory. With data from 3.6 to 500 μm, dust models are strongly constrained. Using a physical dust model, for each pixel in each galaxy we estimate (1) dust surface density, (2) dust mass fraction in polycyclic aromatic hydrocarbons (PAHs), (3) distribution of starlight intensities heating the dust, (4) total infrared (IR) luminosity emitted by the dust, and (5) IR luminosity originating in subregions with high starlight intensity. The dust models successfully reproduce the observed global and resolved spectral energy distributions. With the angular resolution of Herschel, we obtain well-resolved maps (available online) for the dust properties. As in previous studies, we find the PAH fraction q PAH to be an increasing function of metallicity, with a threshold oxygen abundance Z/Z e ≈0.1, but we find the data to be fitted best with q PAH increasing linearly with log O H ( ) above a threshold value of 0.15(O/H) e . We obtain total dust masses for each galaxy by summing the dust mass over the individual map pixels; these "resolved" dust masses are consistent with the masses inferred from a model fit to the global photometry. The global dust-to-gas ratios obtained from this study are found to correlate with galaxy metallicities. Systems with Z/Z e 0.5 have most of their refractory elements locked up in dust, whereas in systems with Z/Z e 0.3 most of these elements tend to remain in the gas phase. Within galaxies, we find that q PAH is suppressed in regions with unusually warm dust with n m n  L L 70 m 0.4 dust ( ) . With knowledge of one long-wavelength flux density ratio (e.g., f 160 /f 500 ), the minimum starlight intensity heating the dust (U min ) can be estimated to within ∼50%, despite a variation in U min of more than two orders of magnitude. For the adopted dust model, dust masses can be estimated to within∼0.2 dex accuracy using the f 160 /f 500 flux ratio and the integrated dust luminosity, and to∼0.07 dex accuracy using the 500 μm luminosity n m n L 500 m ( )alone. There are additional systematic errors arising from the choice of dust

Introduction
Interstellar dust affects the appearance of galaxies by attenuating short-wavelength radiation from stars and ionized gas and contributing infrared (IR), submillimeter, millimeter, and microwave emission. Dust is also an important agent in the fluid dynamics, chemistry, heating, cooling, and even ionization balance in some interstellar regions, with a major role in the process of star formation. Despite the importance of dust, determination of the physical properties of interstellar dust grains has been a challenging task (for a review, see Draine 2003). Even the overall amount of dust present in other galaxies has often been very uncertain.
The "Key Insights on Nearby Galaxies: A Far-Infrared Survey with Herschel" (KINGFISH; Kennicutt et al. 2011) project is an imaging and spectroscopic survey of 61 nearby (distance < D 30 Mpc) galaxies with the Herschel Space Observatory. The KINGFISH galaxy sample was chosen to cover a wide range of integrated properties and local interstellar medium (ISM) environments found in the nearby universe. KINGFISH is a direct descendant of the "Spitzer Infrared Nearby Galaxies Survey" (SINGS; Kennicutt et al. 2003), which produced complete Spitzer imaging with the Infrared Array Camera (IRAC; Fazio et al. 2004) and the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) instruments on Spitzer Space Telescope (Werner et al. 2004). The new Herschel observations include a complete mapping of the galaxies with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Spectral and Photometric Imaging Receiver (SPIRE; Griffin et al. 2010) instruments. The merged KINGFISH and SINGS data set provides panchromatic mapping of the galaxies, across a wide range of local extragalactic ISM environments. In addition, we have KINGFISH and SINGS data for nine additional galaxies that fell within the 61 KINGFISH target fields. The photometric maps cover wavelengths from 3.6 to 500 μm, allowing us to produce well-resolved maps of the dust in nearby galaxies. Skibba et al. (2011) modeled the dust in the KINGFISH galaxy sample using "modified blackbody" models. In the present work, we employ a physically motivated dust model based on a mixture of amorphous silicate grains and carbonaceous grains, each with a distribution of grain sizes (Draine & Li 2007, hereafter DL07). The dust grains are heated by starlight, and the model allows for a distribution of intensities for the starlight heating the dust. With a small number of adjustable parameters, the DL07 model reproduces the observed spectral energy distribution (SED) of the dust emission for a variety of astrophysical systems, giving some confidence in the reliability of dust masses estimated using the model. The DL07 model has been found to be consistent with the 3.6-500 μm emission from dust in the star-forming galaxies NGC628 and NGC6946 , dust across M31 (Draine et al. 2014), emission from annular rings in the KINGFISH galaxy sample (Hunt et al. 2015), and overall dust SEDs from KINGFISH galaxies (Dale et al. 2017).
The present work is a sequel to the KINGFISH study of NGC628 and NGC6946 (Aniano et al. 2012, hereafter AD12). AD12 developed the image processing and dust modeling techniques employed here, using the spiral galaxies NGC628 and NGC6946 as examples. The present work takes into account a recent "recalibration" of the DL07 model made possible by Planck observations of diffuse Galactic emission (Planck Collaboration et al. 2016). We expand the spatially resolved dust modeling to the full KINGFISH galaxy sample, producing maps of dust mass surface density, polycyclic aromatic hydrocarbon (PAH) fraction, and intensities of the starlight heating the dust. Dependences of dust/ gas ratio and PAH abundance on galaxy metallicity are examined, and resolved trends within galaxies are studied. While the present results are undoubtedly model-dependent, comparison of different dust models is beyond the scope of the present work.
The paper is organized as follows. A brief overview of the KINGFISH sample is given in Section 2, and in Section 3 we discuss the data sources. Background subtraction and data processing are described in Section 4, and the dust model is summarized in Section 5, including the Planck-based dust mass "recalibration" (Section 5.2). Results are reported in Section 6 with a comparison of dust parameter estimates based on different dust modeling strategies given in Section 6.4; global trends with metallicity are described in Sections 6.5 and 6.6; and resolved trends of DL07 fitted parameters are discussed in Section 6.7. We summarize the main results in Section 7. Appendix A (online version) displays maps of selected dust parameters for each of the 62 galaxies where we have reliable dust detections, at both MIPS160 and SPIRE250 resolution. Appendix B describes the method used to obtain upper limits for the dust mass for the eight galaxies (five dwarfs, three ellipticals) where we were unable to measure the dust mass reliably. The online data set with the KINGFISH data and dust models is described in Appendix C. In Appendix D we examine the robustness of the results as the point-spread function (PSF) is reduced, precluding use of the lower-resolution cameras (e.g., MIPS160 and SPIRE500).

Galaxy Sample
The observational program was designed to cover the 61 galaxies in the KINGFISH galaxy sample. Because we will also be discussing the nine extra galaxies and the statistical properties of various subsamples, we list these for clarity in Table 1. For each galaxy, we list in Table 2 the type, adopted distance, and major and minor optical radii (corresponding to ∼25th mag arcsec −2 isophotes), all taken from Kennicutt et al. (2011, Table 1).
The galaxies IC 3583, NGC 586, NGC 1317, NGC 1481, NGC 1510 were not part of the KINGFISH sample, but were observed because each happened to be in the field of view of a KINGFISH galaxy. For these galaxies, we have our standard imaging with PACS and SPIRE, as well as prior observations with IRAC and MIPS, so we are able to measure and model their SEDs with the same techniques as the KINGFISH galaxies. Information for these nine "extra" galaxies is appended to many of the tables below.

Metallicities
Table 2 also lists the oxygen abundance 12+log 10 (O/H) for the galaxies in our sample. These are "characteristic" abundances, which Moustakas et al. (2010) take to be the values at galactocentric radius R=0.4R 25 . For six of the KINGFISH galaxies (DDO 154, IC 342, NGC 628, NGC 2146, NGC 3077, and NGC 5457), we use metallicities based on observations of weak lines (specifically, [N II]5726 and [O III]4364) that allow "direct" determination of the electron temperature in the H II regions responsible for the line emission (Storchi-Bergmann et al. 1994;van Zee et al. 1997;Pilyugin et al. 2007;Engelbracht et al. 2008;Li et al. 2013;Croxall et al. 2016). For these galaxies, we list the preferred weak-line metallicities in the PP04N2 column.
For the remaining 55 KINGFISH galaxies, we consider two popular "strong line" estimators: the "PT" or Pilyugin & Thuan (2005) method, taken from Moustakas et al. (2010), and the "PP04N2" method, based on [NII]/Hα (Pettini & Pagel 2004). Abundance measurements by Moustakas et al. (2010, "char-acteristic" values from their Table 8) with the "KK04" (Kobulnicky & Kewley 2004) calibration were converted to PP04N2 values, according to the parameters recommended by Kewley & Ellison (2008). This procedure is described in detail by Hunt et al. (2016), who use the same metallicities in their analysis; they preferred the PP04N2 calibration because it shows tighter scaling relations overall than other calibrations, and because its behavior in the mass-metallicity relation is quite similar to weak-line electron temperature determinations (e.g., Andrews & Martini 2013).
For NGC 1316, NGC 2841, and NGC 5055, the original KK04 O/H values (∼9.4) from Moustakas et al. (2010) exceeded the range of applicability for the transformations formulated by Kewley & Ellison (2008). Thus we have (somewhat arbitrarily) given these three galaxies a maximum metallicity of 12+log 10 (O/H)=9.0, consistent with what is advocated by Pilyugin et al. (2007). Ultimately, the metallicities for these three galaxies are uncertain, but toward the high end of the observed range. Figure 1 compares the PT and PP04N2 metallicity estimates for the 55 KINGFISH galaxies where "direct" method estimates are unavailable. Note that the PT and PP04N2 metallicities differ by as much as 0.5 dex (e.g., DDO 154, type IBm) or even 0.63 dex (e.g., NGC 1482, type SA0). It is evident that the metallicity estimates have significant uncertainties, and that there are systematic differences between the two methods (see also Kewley & Ellison 2008). Below we will argue, by comparing PAH abundances estimated from infrared observations with these two metallicity estimates, that the PP04N2 estimate appears to be more reliable, at least for the galaxies in the KINGFISH sample.

Infrared, Far-infrared, and Submillimeter
Most of the galaxies in the KINGFISH sample are part of the SINGS galaxy sample and were imaged by Spitzer Space Telescope as part of the SINGS observing program (Kennicutt et al. 2003). IRAC and MIPS imaging obtained by other Spitzer Space Telescope observing programs was available for the remaining KINGFISH galaxies.
The KINGFISH project imaged the galaxies with the Herschel Space Observatory (Pilbratt et al. 2010), following the observing strategy described by Kennicutt et al. (2011), using the 70, 100, and 160 μm PACS filters and the 250, 350, and 500 μm SPIRE filters. The maps were designed to cover a region out to 1.5 times the optical radius R 25 , with good signal-to-noise ratio (S/N) and redundancy.
Following AD12, we will use "camera" to identify each optical configuration of the observing instruments; that is, each different channel or filter arrangement of the instruments will be referred to as a different "camera." With this nomenclature, each "camera" has a characteristic spectral response and PSF. We will refer to the IRAC, MIPS, PACS, and SPIRE cameras using their nominal wavelengths in microns: IRAC3.6, IRAC4.5, IRAC5.8, IRAC8.0, MIPS24, MIPS70, MIPS160, PACS70, PACS100, PACS160, SPIRE250, SPIRE350, and SPIRE500.
IRAC imaged the galaxies in four bands, centered at 3.6, 4.5, 5.8, and 8.0 μm, as described by Kennicutt et al. (2003). The images were processed by the SINGS Fifth Data Delivery pipeline. 32 The IRAC images are calibrated for point sources. Photometry of extended sources requires so-called "aperture corrections." We multiply the intensities in each pixel by the asymptotic (infinite radii) value of the aperture correction (i.e., the aperture correction corresponding to an infinite radius aperture). We use the factors 0.91, 0.94, 0.66, and 0.74 for the 3.6 μm, 4.5 μm, 5.8 μm, and 8.0 μm bands, respectively, as described in the IRAC Instrument Handbook (V2.0.1). 33 Imaging with MIPS at 24, 70, and m 160 m was carried out following the observing strategy described in Kennicutt et al. (2003). The data were reduced using the Local Volume Legacy (LVL) project pipeline. 34 A correction for nonlinearities in the MIPS70 camera was applied, as described by Dale et al. (2009) and Gordon et al. (2011).
The galaxies were observed with the PACS and SPIRE instruments on Herschel, using the "Scan Map" observing mode. Both PACS and SPIRE images were first reduced to "level 1" (flux-calibrated brightness time series, with attached sky coordinates) using HIPE v11.1.0 (Ott 2010), and maps ("level 2") were created using the Scanamorphos data reduction pipeline (Roussel 2013), v24.0. This reduction strategy used  the latest available PACS and SPIRE calibrations (as of 2014 July) and was designed to preserve the low-surface-brightness diffuse emission. The assumed beam sizes are 465.4, 822.6, and 1769 arcsec 2 for SPIRE250, SPIRE350, and SPIRE500, respectively. Additionally, we excluded discrepant bolometers from the map and adjusted the pointing to match the MIPS24 map.

H I Observations
To measure the H I gas mass, we use H I 21 cm line observations made with the NSF's NRAO 35 Karl G. Jansky Very Large Array (VLA).
For 23 of our galaxies, we have data from The H I Nearby Galaxies Survey (THINGS; Walter et al. 2008), and for four galaxies we use data from the LittleTHINGS survey (Hunter et al. 2012). For 10 galaxies without THINGS or Little-THINGS observations, we obtained VLA 21 cm maps in programs AL731 and AL735, in some cases also incorporating archival VLA observations. For eight targets, we reduced and incorporated VLA archival observations of the 21 cm line. For one galaxy, NGC4559, we use archival Westerbork Synthesis Radio Telescope observations. These observations are described in Leroy et al. (2013). For each galaxy, the source of the H I map is listed in Table 3. The dominant uncertainty on the measured H I masses comes from the calibration uncertainties of ∼10%.
For NGC 1266, Alatalo et al. (2011) 10 6  based on 21 cm absorption of the radio continuum from the nucleus (for an assumed T spin = 100 K). However, we estimate that H I 21 cm line emission from as much as ∼2×10 9 M  of H I could have gone undetected because of the strong continuum (0.1 Jy at 1.4 GHz), so the H I mass must be considered highly uncertain.
Thus we have H I data for 57 of the 61 KINGFISH galaxies. H I 21 cm observations were not available for NGC 1316 (SAB0), NGC 1377 (S0), NGC 1404 (E1), or NGC 5866 (S0), nor for the nine extra galaxies.

CO Observations
To estimate H 2 masses, we use observations of CO line emission together with an assumed ratio of H 2 mass to CO luminosity. The adopted CO-to-H 2 "conversion factors" are discussed in Section 4.5.
For 38 KINGFISH galaxies, we use 12 CO = J 2 1 maps from the the HERA CO Line Emission Survey (HERACLES; Leroy et al. 2009Leroy et al. , 2013. For NGC4826, we use 12 CO = J 1 0 mapping from the Nobeyama Radio Observatory (J. Koda et al. 2019, in preparation). We propagate uncertainties on the CO integrated intensities from the spectra through the gridding and masking of the cube as described in Leroy et al. (2013).
Thus we have CO data for 41 of the 61 galaxies in the KINGFISH sample. CO observations are not available for any of the nine extra galaxies.

Background Subtraction
All camera images are first rotated to RA/Dec coordinates and then trimmed to a common sky region. For each image, we estimate the best-fit "tilted plane" background (consisting of instrumental background, Galactic foreground emission, and cosmic infrared background emission) using an iterative procedure described in AD12. The procedure uses multiple cameras to identify regions in the image where only background emission is present. Regions where excess emission is detected at more than one wavelength are not used for background estimation.

Convolution to Common Resolution
After background subtraction, the images are convolved to a common PSF and resampled on a common final-map grid, with pixel sizes for each final-map PSF as given in Table 4. Finally, the dispersion in intensities of the background pixels (which includes noise coming from unresolved, undetected background sources) is used to estimate the pixel flux uncertainties. By comparing the MIPS and PACS images, we can also estimate a calibration uncertainty. The procedures used are fully described in AD12.
As discussed in AD12, multiwavelength observations must be degraded to a common PSF before dust models are fit to the observed intensities. The convolution to a common PSF is carried out using the methods described by Aniano et al. (2011). In the present work, we present, for each galaxy, resolved results at two final-map PSFs: SPIRE250 and MIPS160, henceforth abbreviated as S250 and M160. S250 is the PSF with the smallest FWHM that allows use of enough cameras to adequately constrain the dust SED (IRAC, MIPS24, PACS70, 100, 160, and SPIRE250). The M160 PSF allows inclusion of all the cameras (IRAC, MIPS, PACS, SPIRE), therefore producing the most reliable maps; this will be our "gold standard." Table 4 lists the resolutions of the cameras, the pixel size in the final-map grids used, and the other cameras  that can be used at this resolution. In Appendix D we compare dust mass estimates obtained with different final-map PSFs.

Image Segmentation
After convolution to a common "final-map" PSF and background subtraction, we next fit a dust model to the observed SED of each pixel in the field. In order for dust mass estimation to be reliable, the pixel's SED must be measured in a number of bands with a reasonable S/N. However, estimates of the total dust infrared luminosity per unit area from a single pixel, S Ld , are reliable so long as there is a significant detection of far-infrared emission after background subtraction.
The procedure used for automatically identifying "galaxy" pixels is described in Appendices A and B of AD12. For purposes of dust mass estimation, we need to limit the modeling to a "galaxy mask" consisting of pixels where the emission from the galaxy of interest has sufficiently high surface brightness for dust mass estimation (via SED fitting) to be reasonably reliable.
A simple criterion for "sufficiently high surface brightness" is that the total dust luminosity/projected area S Ld exceed a specified threshold value, S Ld,min . The value chosen for S Ld,min will depend on the noisiness of the data (which may depend on the brightness of the (subtracted) Galactic foreground emission, as well as on the presence of other extragalactic objects in the field, stars, or even small-scale structure in the Galactic foreground, which may compromise background estimation and subtraction). The choice of S Ld,min will also depend on the choice of final-map PSF: use of a larger PSF improves the S/N in each pixel by smoothing and also enables use of more cameras to constrain the dust modeling, and thus may allow use of a lower threshold S Ld,min . In the present study, S Ld,min was chosen subjectively for each galaxy.
In this paper, we report results for two final-map PSFs, S250 and M160, with the PSF FWHM corresponding to linear scales = FWHM 0.88 kpc for S250 and1.88 kpc for M160 at the median distance = D 10 Mpc of the KINGFISH galaxy sample. Table 2 lists S Ld,min used to define the galaxy masks for the M160 resolution studies, for the 62 galaxies where we detect dust emission. Our adopted values of S Ld,min vary from galaxy to galaxy, ranging from values as low as For each galaxy where dust is reliably detected at M160 resolution, we also generate a S250-resolution mask, intended to comprise the region where the S250-resolution data permit reliable estimation of the dust surface density. Our S250 masks are often similar in size to the M160 mask, but for some galaxies the S250 mask is considerably smaller than the M160 mask; the most extreme example is NGC 1481, where the S250 mask area is only 37% of the M160 mask area. The M160 and S250 masks are shown in . The solid angle of each mask is listed in Table 2. Because of the improved S/N in each pixel, most of the analysis in this paper will be done with the M160-resolution images and masks.
For five dwarf galaxies where dust detection is uncertain (DDO 053, DDO 154, DDO 165, Hol1, and M81dwB), we choose instead to use masks defined by H I observations. For three elliptical galaxies where dust detection is uncertain (NGC 0584, NGC 0855, and NGC 1404), we useS Ld,min -based masks.We do not detect dust in any of these eight galaxies. See Appendix B for further details.

Integrated Fluxes
The Spitzer and Herschel band surface brightnesses are integrated over the M160-and S250-resolution galaxy masks to obtain integrated flux densities. The IRAC and MIPS flux densities are given in Table 5, and the PACS and SPIRE flux densities are given in Table 6. Note that MIPS70, MIPS160, SPIRE350, and SPIRE500 are not used at S250 resolution.
The uncertainties given in Tables 5 and 6 include uncertainties associated with background subtraction, as well as calibration uncertainties. As discussed in Section 5.3, the fluxes measured by PACS and MIPS sometimes differ by considerably more than the estimated uncertainties: we know that some of the uncertainties have been underestimated, although it is not clear how to improve on our estimates. Dale et al. (2012Dale et al. ( , 2017 carried out careful foreground star and background galaxy removal tailored for globally integrated photometry. For 44 of the 53 KINGFISH galaxies where we

10
The Astrophysical Journal, 889:150 (39pp), 2020 February 1 Aniano et al. claim dust detections, the SPIRE500 flux for our M160 galaxy mask is within 10% of the global SPIRE500 photometry from Dale et al. (2017). Thus, we are not missing a significant reservoir of dust in the outer parts of the disk. 36

Gas Masses
For galaxies observed by THINGS, H I 21 cm line intensities were extracted over the area of the M160-resolution galaxy mask for each galaxy. The H I column density N H I ( ) was estimated assuming the 21 cm emission to be optically thin.
For the 38 galaxies in the HERACLES sample, 12 CO(2-1) line fluxes were obtained by integrating over the M160-resolution galaxy mask, and the H 2 mass was estimated from the CO  2 1 line flux assuming and a standard conversion factor 37 =----X 2 10 H cm K km s CO,1 0 20 2 2 1 1 ( ) . The adopted X CO value is representative of the values found in 26 nearby star-forming galaxies by Sandstrom et al. (2013).
For NGC 1266, we use the integrated CO emission and the lower bound on the H I mass from Alatalo et al. (2011).

DL07 Dust Model
We employ the DL07 dust model, using "Milky Way" grain size distributions (Weingartner & Draine 2001a). DL07 described the construction of the dust model, and AD12 described its usage in the context of the KINGFISH galaxies. The DL07 dust model has a mixture of amorphous silicate grains and carbonaceous grains, with a distribution of grain sizes. The distribution of grain sizes was chosen to reproduce the wavelength dependence of interstellar extinction within a few kiloparsecs of the Sun (Weingartner & Draine 2001a). The silicate and carbonaceous content of the dust grains was constrained by observations of the gas phase depletions in the ISM. It is assumed that the radiation field heating the dust has a universal spectrum, taken to be that of the local interstellar radiation field estimated by Mathis et al. (1983), scaled by a dimensionless factor U. Following DL07, we assume that in each pixel there is dust exposed to radiation with a single intensity U min , and also dust heated by a power-law distribution of starlight intensities with < < U U U min max : where M d is the total dust mass in the pixel, and γ is the fraction of the dust mass that is heated by the power-law distribution of starlight intensities. The DL07 model has six adjustable parameters pertaining to the dust and the starlight heating the dust: 1. q PAH : fraction of the total grain mass contributed by PAHs containing fewer than 10 3 carbon atoms. 2. U min : intensity of the diffuse ISM radiation field heating the dust, relative to the solar neighborhood. 3. α: exponent of the power-law distribution of heating starlight intensities between U min and U max . The case α=2 corresponds to constant dust heating power per logarithmic interval in starlight intensity U; many galaxies seem to be characterized by α≈2. 4. U max : maximum heating starlight intensity of the powerlaw distribution of heating starlight intensities. 5. γ: fraction of the dust mass exposed to the power-law distribution of starlight intensities. 6. M d : dust mass in the pixel.
In addition, for modeling the observed fluxes in the IRAC bands, we have an additional adjustable parameter (see AD12): 7. Ω å : solid angle subtended by stars within the pixel, determined from the "direct" starlight intensity in the infrared, that is, starlight that directly contributes to the IRAC photometry, without warming the dust. The stars are assumed to have color temperature 5000 K.
The mean starlight intensity seen by the dust is The parameter γ is directly related to f PDR , defined to be the fraction of the total dust luminosity L d that is radiated by dust in regions where > U 10 2 , typically photodissociation regions (PDRS) near luminous stars: max ( ) and the adopted grain size distribution and grain properties, the dust emission spectrum is computed from first principles. The observed SEDs are consistent with models having = U 10 max 7 , and we therefore fix º U 10 max 7 . Moreover, the model emission is linear in M d ,  L , and γ (or, equivalently, f PDR ), so in the dust fitting algorithms we only need to explore a three-dimensional parameter space (q PAH , U min , and α). The limits on adjustable parameters are given in Table 7. The allowed range for U min is determined by the wavelength coverage of the data used in the fit.
The region observed is at a distance D from the observer, and Ω j is the solid angle of pixel j. For each pixel j, the best-fit model vector ). For each pixel j, we find the best-fit model parameters { } by minimizing c 2 , as described by AD12. After the resolved (pixel-by-pixel) modeling of the galaxy is performed, we compute a set of global quantities by adding or taking weighted means (denoted as á ñ ... ) of the quantities in each individual pixel of the map. The total dust mass M d , total dust luminosity L d , and total dust luminosity radiated by dust in regions with > U 10 2 , L d,tot , are given by where the sums extend over all the pixels j that correspond to the target galaxy (i.e., the "galaxy mask" pixels, as described in AD12). The dust-mass-weighted PAH mass fraction á ñ q PAH and mean starlight intensity á ñ U are given by

( )
We similarly define the dust-mass-weighted minimum starlight intensity While the average value of α is of little physical significance (the sum of two power laws is not a power law), for the purpose of discussion we define a representative value

( )
We also fit a dust model to the global photometry of each galaxy (i.e., a single-pixel dust model). Below we will compare the result of this single-pixel global model with summing over the fits to individual pixels.

Post-Planck Renormalization of DL07 Dust Masses and Starlight Intensities
Planck Collaboration et al. (2016) fitted the DL07 dust model to all-sky maps in the Planck 857, 545, 353, 217, 143, and 100 GHz (350μm, 550μm, 850μm, 1.4 mm, 2.1 mm, and 3.0 mm) bands; DIRBE 100μm, 140μm, and 240μm bands; IRAS 60 and 100μm bands; and the WISE m 12 m band, to estimate the dust mass surface density for over 50 million ¢´¢ 1.7 1.7 pixels. About 270,000 of these pixels contain spectroscopically confirmed Sloan Digital Sky Survey (SDSS) quasars, which were used to estimate the correlation of quasar reddening with the reddening predicted by the DL07 dust model. It was discovered that the DL07 model tends to overpredict the reddening by a factor of ∼2. The Panchromatic Hubble Andromeda Treasury study of stars in M31 (Dalcanton et al. 2015) also found that the DL07 dust model, if constrained to reproduce the observed infrared emission (Draine et al. 2014), overpredicted the reddening of stars in M31 by a factor of ∼2. The SDSS quasars allow the bias factor to be estimated: Planck Collaboration et al. (2016) found that the bias appeared to depend on the value of U min : If the reddening -E B V ( ) has been overestimated, it is reasonable to suppose that the dust mass per area has also been overestimated, by approximately the same factor as the reddening. Therefore, we will correct the DL07 dust mass estimates by the same empirical correction factor as for the reddening.  Because the KINGFISH sample includes some pixels with high U min , we choose to limit the Planck-derived correction factor for > U 1 min : Because the dust models are required to reproduce the dust luminosity, a reduction in the estimated amount of dust implies a corresponding increase in the estimated starlight intensities. Thus we take, for pixel j, All of the dust and starlight parameters (M d , U min ) reported below are "renormalized" values from Equations (13) and (16). The above "renormalization" is required because observations indicate that the far-infrared and submillimeter opacity of interstellar dust per unit reddening is somewhat larger than the DL07 model values. 38 We will find below ( Table 9) that the global average á ñ C dust for the 62 galaxies ranges from 0.45 to 0.69, with median á ñ = C 0.62 dust . Thus the typical correction relative to the DL07 model is a reduction in M d by a factor of ∼0.62.
Because Equation ( , the estimated dust-to-gas ratios (DGR) in regions with  U 1 min are reduced by a constant factor of 0.70. This applies to the study of Sandstrom et al. (2013), which was dominated by regions with á ñ > U 1. However, because = C constant dust , the CO-to-H 2 ratios found by Sandstrom et al. (2013) are unaffected by the renormalization.

Why Both MIPS and PACS Are Needed
As discussed by Aniano et al. (2011), the MIPS160 PSF cannot be convolved safely into any of the PSFs of the remaining cameras. Therefore, if we wish to include MIPS160 photometry in the dust modeling, we must "degrade" all other images into the MIPS160 PSF.
There are two reasons why we want to include MIPS160 even though PACS160 imaging is available. First, using the larger PSF increases the S/N for the imaging, thereby allowing photometry to be extended to lower surface brightness regions. Second, there are significant and unexplained discrepancies between PACS160 and MIPS160 photometry. Similar discrepancies are found between PACS70 and MIPS70. Figure 2 shows histograms of the global PACS70/MIPS70 flux ratio (left panel) and the global PACS160/MIPS160 flux ratio (right panel) for each of the KF62 galaxies with reliable dust detections. Each histogram shows the names of galaxies in the bin: "NGC," "DDO," "Holmberg," and "IC" are abbreviated to "n," "d," "Hol," and "i," respectively.
The PACS70 and MIPS70 bandpasses differ slightly, as do the PACS160 and MIPS160 bandpasses. However, AD12 show that for reasonable dust SEDs, the slight difference in bandpasses can explain differences in reported fluxes of only 9% at 70 μm and ∼2% at 160 μm, whereas much larger PACS/MIPS discrepancies are often observed.
AD12 ( Table 1). Eight galaxies where dust was not reliably detected have been excluded (see text).  Figure 2 illustrates that even after summing over the full galaxy mask, PACS70 and MIPS70 often disagree by more than a factor of 1.2, and sometimes up to a factor of 1.4. The median ratio is 1.17.
PACS160 and MIPS160 are generally in better agreement, but often have discrepancies larger than 10%. There are two outliers in Figure 2: NGC 2146 (PACS160/MIPS160=1.6) and NGC 2915 (PACS160/MIPS160=2.3). The high value of PACS160/MIPS160 for NGC 2146 may be the result of the sublinear response of MIPS160 on the very bright nucleus of NGC 2146. The case of NGC 2915 is unclear: the peak surface brightness is modest. Perhaps the background has been oversubtracted in the MIPS160 image or undersubtracted in the PACS160 image.
Because it is usually unclear why PACS and MIPS disagree (the discrepancies are too large to be attributed to differences in bandpasses), we consider that both PACS and MIPS photometry should be included if we wish to estimate the dust parameters with the best accuracy available. AD12 also found that, for a given camera set, dust parameter estimates do not change significantly when using a broader PSF, so modeling at MIPS160 PSF does not significantly alter the dust parameter estimates. We consider our "gold standard" (i.e., the PSF and camera combination that gives the most accurate dust parameter estimates) to be resolved (i.e., multipixel) modeling done using the MIPS160 PSF, using photometry from all of the IRAC, MIPS, PACS, and SPIRE cameras.

Results
For each galaxy in the KF62 sample, Table 9 presents the global dust parameters estimated for the "gold standard" modeling, including information characterizing the intensity of the starlight heating the dust in each galaxy. The modeling was done at MIPS160 PSF, using all the cameras available; we also give results of modeling at S250 resolution.
The given quantities are obtained by summing or averaging over the resolved maps using Equations (7)-(10). The dust masses listed in Table 9 were obtained using the DL07 model, but then "renormalized" following Equation (13). The renormalization factor C d depends on U min and therefore varies from pixel to pixel. The overall renormalization factor for each galaxy is given in Table 9, for both M160 and S250 resolution. Henceforth, M d , U min , and U will refer to the renormalized values of these quantities (see Equations (13)- (16)).

One Example: NGC 5457=M101
To illustrate the quality of the data and the modeling results for the KINGFISH galaxies, we choose the large, nearly faceon spiral NGC 5457 (M101) as an example. As for all our galaxies, the dust mass, PAH abundance, and starlight heating parameters are adjusted separately for each pixel.
The parameter α characterizes the distribution of starlight intensities heating dust within a pixel (see Equation (1)). Figure 3 shows maps of the best-fit α values for the M160-and S250-resolution modeling. At M160 resolution, α is azimuthally coherent but has a notable radial gradient, with a » 1.7 in the center, and α≈2.3 beyond galactocentric radius ∼6 kpc.
While the variation in best-fit α is apparent, these values are all close to α=2, the case where there is equal power per unit U log . At S250 resolution, the S/Ns are lower, and the S350, S500, and M160 cameras are not used; the α map for the S250resolution modeling shows more pixel-scale variations, but with a radial trend similar to the M160-resolution modeling.
In general, the DL07 model successfully reproduces the resolved SEDs in M101. Figure 4 compares the model 500 μm surface brightness with observations. The upper panel shows modeling at M160 resolution (the observed SPIRE500 intensity is used as a model constraint). The DL07 model is generally within ±10% of the observed SPIRE500 intensity, except at the outer edges of the mask where the S/N is low. The model appears to fall short by ∼10% in the outer regions (galactocentric radius ∼15 kpc=0°.13), where the metallicity has dropped to + » 12 log O H 8.25 10 ( ) (Li et al. 2013). This could indicate that the frequency dependence of the dust opacity becomes less steep as the metallicity drops, which is consistent with the SED of the SMC (Bot et al. 2010;Israel et al. 2010;Planck Collaboration et al. 2011;Draine & Hensley 2012) and with evidence for a submillimeter excess in galaxies with metallicities The lower panel of Figure 4 compares modeling at S250 resolution (no data longward of 250 μm used to constrain the model) with the SPIRE500 observations. In the bright spiral arms, the 500 μm intensity is overpredicted by ∼25%. Once again we see a radial gradient: the model overpredicts SPIRE500 in the central regions and underpredicts SPIRE500 at =   R 8 kpc 0 . 07. In the outer regions the fit is poorer, presumably due to the low S/N at S250 resolution. Figure 5 shows maps of dust and starlight heating parameters for M101. There are two sets of figures: the first set (rows 1 and 2) corresponds to modeling done at M160 resolution, using data from all (IRAC, MIPS, PACS, and SPIRE) cameras, that is, "gold standard" modeling, and the second (rows 3 and 4) corresponds to modeling done at S250 resolution, using IRAC, MIPS24, PACS, and SPIRE250 cameras. This latter modeling is able to resolve smaller scale structures in the galaxies, but is overall less reliable, particularly in the outer regions where the surface brightness is lower and dust is cooler.
Because of the proximity of M101 ( = D 6.7 Mpc), the spiral structure is visible even at M160 resolution. At M160 resolution At S250 resolution, the peak at NGC 5461 has a dust/ luminosity/area S = -L 10 pc Ld 3.2 2  (corresponding to a dust luminosity =Ĺ L 6 10 d 7  in a single 195×195 pc 2 S250 map pixel). Thus at S250 resolution, we are able to measure the IR emission from the dust over a dynamic range of ∼2000 in S Ld .
Maps of dust surface density S Md are also shown for both the M160 and S250 modeling. At both M160 and S250 resolution, S Md has a peak at the extranuclear luminosity peak. At S250 resolution, we estimate a peak dust surface density5 Maps of the starlight modeling parameter U min,DL07 are also shown at both M160 and S250 resolution. In M101, U min,DL07 ranges from values as high as 30 (the largest value permitted by our modeling) to values as low as ∼0.07 in the outer parts of the galaxy. The highest values of = U 30 min arise in the S250 modeling, with high values of U min appearing in a fraction of pixels in low-surface-brightness regions to the east of the center. The high U min values found in these regions using S250resolution data are probably unphysical, arising as the result of low-S/N data: an upward fluctuation in PACS70 (or a downward fluctuation in SPIRE250) can drive the fitting to a high U min value. Within ∼5 kpc of the center, with higher surface brightnesses, we generally find   U 0.5 4 min . And in the M160 modeling, we do not obtain very high values of U min even in the low-surface-brightness outer regions.
Maps of q PAH are also shown at both M160 and S250 resolution. The modeling finds a very high value of q PAH along the SSE edge of the galaxy; this is seen in both the M160 and S250 modeling of an extended region approximately 12 kpc SSE of the center. The high estimates for q PAH could arise from errors in the IRAC 5.6 and 8μm photometry, probably due to errors in background subtraction. . At M160 resolution, where all cameras are used to constrain the model, α is azimuthally coherent but with a radial gradient: a » 1.7 in the center, and a » 2.3 in the outer regions. The S250 map is noisier, because not all of the cameras can be used, and the S/N of the bands that can be used is reduced.  Maps of f PDR -the fraction of the dust luminosity that is contributed by dust heated by starlight with > U 100-are shown at both M160 and S250 resolution. High values of f PDR are seen at many of the positions where S Ld peaks, which is consistent with the idea that these are regions with active star formation, with some fraction of the dust exposed to intense radiation fields in or near OB associations. However, we also see high f PDR values in some of the lowest-surface-brightness regions near the edge of the galaxy mask; this is presumably an indication that photometric errors and errors in background subtraction are leading to overestimation of 24 or m 70 m emission relative to the total dust luminosity. Thus our derived values of f PDR appear to be unreliable in the lowest-surfacebrightness regions.
We also show the global SED for M101, extracted from the galaxy mask. In the upper right panel, the rectangular symbols show the measured fluxes ±1σ for the seven Spitzer cameras and the six Herschel cameras. At 70 and 160μm, both red and blue rectangles are shown, with the MIPS and PACS photometry. Also shown is a single-pixel DL07 model, where the DL07 model is fitted to the global photometry. The diamonds show the model fluxes for each of the instrumental bandpasses. In the case of M101, the model (with six adjustable parameters: ) is consistent with the photometry at 11 independent wavelengths (3.6μm to 500μm). In row 3, column 3, we show a single-pixel DL07 model fitted to only the photometry that is used for the S250 modeling (i.e., MIPS70, MIPS160, SPIRE350, or SPIRE500 are not used when adjusting the model parameters). The dashed rectangles show these unused measurements; we see that for M101 the single-pixel model does quite well at predicting the fluxes at 350 and 500μm, with only a 1σ underprediction even at 500μm. The single-pixel global fit parameters are given in the SED plots. Table 8 compares total dust mass estimates for M101. Column 2 reports the dust mass estimated from the DL07 model at either M160 or S250 resolution, after summing the dust model over the galaxy mask. Because we opted to use the same S Ld,min for the S250 and M160 modeling, the galaxy masks for the two cases are essentially the same. Column 3 reports the result of fitting a DL07 model to the global photometry; this is referred to as "single pixel" modeling. In columns 4 and 5, we show the multipixel or single-pixel dust masses after renormalizing following Equation (13).
Multipixel versus single-pixel modeling is of course expected to produce different estimates because the models are nonlinear. One notes in Table 8 that the discrepancies between the multipixel and single-pixel mass estimators are reduced when going from the original DL07 model to the renormalized model. It is not clear why this is the case, but this is a welcome result.

Full KINGFISH Sample
Dust is detected reliably for every galaxy in the KF62 sample. Selected images for each of these galaxies are given in Appendix A (Figures 17.1-17.62), following the scheme used for M101 in Figure 5. This is only a fraction of the maps and images that are available online; see Appendix C for a description of the data set.
The "galaxy mask" for each galaxy is shown for both the M160 PSF and the S250 PSF. As for M101, we have opted to use the same S Ld,min for both the M160 and S250 modeling, so the M160 and S250 S Ld,min -based galaxy masks are nearly identical for each galaxy, except for the eight where dust emission is so weak that we treat them as nondetections (see Appendix B). The flux densities F ν measured by Spitzer and Herschel within the M160 and S250 galaxy masks for each galaxy have been given in Tables 5 and 6. The model-derived parameters for the dust and starlight are given in Table 9. Note that α is not included in Table 9 because there is no natural way to define a "mean" α for multipixel modeling. The uncertainties listed for the parameters are based on repeating the fitting procedure with the "observed" fluxes obtained by Monte Carlo sampling from Gaussian distributions with means and widths given by the original observed values and uncertainty estimates (see discussion in Appendix E of AD12). Systematic errors associated with the DL07 model itself have not been estimated. Figures 17.1-17.62 have 12 panels in all, with the top six panels showing results of modeling with the M160 PSF, and the lower six panels repeating this for the S250 PSF. For each PSF, the top row shows maps of the dust luminosity surface density S L d (upper left) and modeled dust surface density S M d (upper center), and the model SED (upper right). The lower row shows maps of the starlight intensity parameter U min (left), the PAH abundance parameter q PAH (center), and the PDR fraction f PDR (right), all restricted to the "galaxy mask." In the SED plot, the observed global photometry is represented by rectangular boxes (Spitzer IRAC and MIPS in red; Herschel PACS and SPIRE in blue). The vertical extent of each box shows the ±1σ uncertainty in the photometry for each band. The black line is the total DL07 model spectrum, and its different components are represented by three colors. The cyan line is the stellar contribution, the dark red line is the emission from dust heated by the power-law U distribution, and the dark green line is emission from dust heated by = U U min . The DL07 model used in this SED plot is a single-pixel model, which tries to reproduce the global photometry by treating the entire galaxy as a single pixel. These "single pixel" models generally do a good job at reproducing the global photometry. Multipixel models, where the photometry in every pixel is fit independently, have many more adjustable parameters, and they naturally do an even better job of reproducing the global photometry after summing over all the pixels in the galaxy mask. It is reasonable to presume that   Note. a Renormalized as described in Section 5.2. models that do a better job of reproducing the photometry will also be preferred for dust mass estimation. The dust mass surface densities and dust luminosities per unit area range over three orders of magnitude in our brightest galaxies.  illustrate that the DL07 model does a satisfactory job of fitting the SEDs. Although each pixel is modeled independently of its neighbors, it is noteworthy that the dust parameters are smoothly varying over the confines of the galaxy, except for the low-surface-brightness outer regions at S250 resolution, where the S/N in individual pixels may become low enough that certain dust and starlight parameters, such as S Md and U min , become somewhat noisy.

NGC 1404
The E1 galaxy NGC 1404 is faint at infrared wavelengths, and there is a foreground star ∼180″ to the SSE. However, NGC 1404 is unambiguously detected in IRAC bands 1-4 and by MIPS24. At λ70μm, there appears to be excess emission at the position of NGC 1404, but the surface brightness is low, and the possibility that the emission is from a foreground or background source cannot be excluded (Dale et al. 2012 NGC 1377 has a compact, dusty core, with an extremely high far-infrared surface brightness. The infrared spectrum (Roussel et al. 2006) shows that it is optically thick at l m  20 m. Weak PAH emission is detected, but because of the uncertain infrared extinction, it is not possible to reliably estimate the PAH abundance parameter q PAH . The dust mass estimates should also be regarded as uncertain because of the unusual nature of the ISM in this galaxy. Figure 6 shows the dust parameter distributions for the 61 KINGFISH galaxies plus nine "extras." The dust parameters shown are the result of the "gold standard" modeling, which is multipixel modeling for each galaxy using the M160 PSF and data from all cameras.

Gold Standard DL07 Fit Results for KINGFISH Galaxies
The first row shows the distributions of M d (left column) and L d (right column) for the KF62 galaxies. The second row shows the distributions of á ñ U (left column) and q PAH (right column), and the bottom row shows the distributions of á ñ f PDR (left column) and a á ñ (right column). In these histograms, the dust masses M d and á ñ U are renormalized, as discussed in Section 5.2. Figure 6 illustrates the large region in the model parameter space spanned by the KINGFISH sample, allowing us to probe the dust properties in a variety of ISM conditions. The total dust mass and dust luminosity found in the galaxies span almost four decades:  (Dellenbusch et al. 2008), NGC 1510 hosts a strong central starburst, and NGC 3049 and NGC 5408 are often classified as starburst galaxies. Thus, the high á ñ f PDR values for these galaxies may be indicative of concentrated star formation or nuclear activity.
Finally, the mean power-law exponent a á ñ spans  1.5 a á ñ  2.5. Allowing α to vary does improve the quality of the fit to the observed SED, but in most cases the fit quality does not suffer greatly if α is held fixed at α=2, reducing the number of free parameters. Recall that α=2 corresponds to equal amounts of dust power per logarithmic interval in starlight intensity U.
Compared with the "gold standard," modeling using PSFs smaller than M160, and hence having fewer cameras available, can affect the derived dust and starlight parameters. As the PSF shrinks, data are provided by fewer cameras, the wavelength coverage shrinks as the PSF is reduced below S500, and the photometry becomes noisier because it is being smoothed over smaller PSFs. Above we have compared two cases: modeling with the M160 PSF versus modeling with the S250 PSF, but additional comparisons are made in Appendix D. Here we simply note some trends. In general, M d is fairly robust: the S250 modeling typically overestimates M d by ∼25%, but agrees with the "gold standard" to within a factor of 1.5 for over 75% of the galaxies (see Figure 19). The q PAH estimates are also robust, with typical changes of less than 15%. Longer wavelength coverage (SPIRE350 and SPIRE500) gives more reliable dust estimates. Even comparing resolved and global modeling of dust properties can give different results; although most parameters are consistent to within a few percent, global modeling can underestimate M d by as much as 35%, as for NGC1481 and NGC3077 (see Figure 21). This is because the resolved models can have "cold" regions with low U min values that contribute to dust mass estimates but do not emerge in the global results (e.g., Galliano et al. 2011;Galametz et al. 2012).

Dependence of Global PAH Fraction on Metallicity
Figure 7(a) shows q PAH versus log O H ( ) for 51 galaxies using direct determinations of (O/H) where available (five galaxies), and PT estimates otherwise. Nineteen galaxies have been omitted: eight dust nondetections have been excluded, NGC 1377 (a dense starburst with a core that is optically thick at 8μm; see Section 6.3.2), plus 10 galaxies for which we have no PT estimate for O/H. The oxygen abundance in these galaxies ranges over more than a factor of 10, and q PAH shows a clear tendency to increase with increasing O/H, although there is considerable scatter. The observed behavior can be approximated by a step function, with an abrupt increase in q PAH when + 12 log O H 10 PT ( ) rises above ∼8.0. Alternatively, q PAH can be approximated by a linear dependence on log O H ( ). The best-fit step function and linear function are shown in Figure 7(a), with c 2 per degree of freedom of 8.0 and 6.6, respectively. Figure 7( ( ) >7.94). This fit, with 53 -2=51 degrees of freedom (dof), has c 2 /dof=3.5: the PP04N2 metallicity is evidently a much better predictor of q PAH than is the PT metallicity. 40 This strongly suggests that the PP04N2 metallicities are more tightly related to the properties of the ISM-including metallicity-that regulate the balance between PAH formation and destruction.
The observed tendency for q PAH to increase with increasing metallicity is consistent with many previous studies. The connection between PAH abundance and metallicity was first noted in ground-based spectroscopy by Roche et al. (1991) and further investigated using Infrared Space Observatory data (Boselli et al. 1998;Madden 2000;Sturm et al. 2000). Hunt et al. (2005Hunt et al. ( , 2010 found PAH emission to be weak in lowmetallicity, blue compact dwarf galaxies. Engelbracht et al. (2005) used IRAC and MIPS24 photometry to show that there was an abrupt drop in the 8 μm/24 μm flux ratio when the metallicity dropped below 8.2, interpreting this as being due to a sharp drop in the abundance of PAHs that normally dominate the emission at 8μm.  estimated q PAH for 61 SINGS galaxies, using the DL07 model with IRAC and MIPS photometry, and found a similar result: a sharp increase in q PAH when + 12 log O H 10 PT ( ) rises above ∼8.2. Nevertheless, there are outliers in Figure 7( Figure 17.8) appears to be reliable. Because the optical spectrum of NGC 1266 is AGN-dominated, the metallicity is not based on emission lines and is therefore highly uncertain. Moustakas et al. (2010) estimated the metallicity from an assumed luminositymetallicity relation. The resulting + 12 log O H 10 PP04N2 ( ) = 8.51 is consistent with the stellar mass-metallicity relation (Andrews & Martini 2013). Perhaps the PAH abundance in this galaxy has been suppressed by phenomena associated with the AGN that is driving a molecular outflow characterized by shocked gas (Alatalo et al. 2011(Alatalo et al. , 2015Pellegrini et al. 2013).
The SAB0 galaxy NGC 1316=Fornax A is another outlier. The dust emission is weak relative to the starlight, making the q PAH estimate uncertain. In addition, the starlight heating the dust is likely from an old population, similar to the bulge of M31, and our estimate of q PAH (based on single-photon heating by starlight assumed to have the solar neighborhood spectrum) Figure 7. (a) PAH abundance parameter q PAH versus oxygen abundance (direct or PT) for 51 galaxies (see text). Two q PAH estimators are shown: one is a step function, and the other is linear above a threshold value. Selected galaxies have been labeled. The step function and linear estimators have similar c = dof 7.8 2 and 6.3. (b) Same but for PP04N2 oxygen abundance, now for 53 galaxies (PP04N2 oxygen abundances are available for NGC 1512 and NGC3077). The five galaxies where O/H has been determined by "direct" methods are shown in green. The linear fit of q PAH versus + 12 log O H 10 ( ) (Equation (18) would then be biased low. The estimate for q PAH in the center of M31 increases by almost a factor of two when calculated using the correct starlight spectrum (Draine et al. 2014), and a similar correction might bring q PAH for NGC 1316 closer to the general trend in Figure 7(b). In addition, the high metallicity estimated for NGC 1316 may be influenced by the AGN contribution to the emission line spectrum.
In Figure 7 it is striking that the bulk of the galaxies with + 12 log O H 10 PP04N2 ( ) >8.3 have q PAH in the 1.5%-5% range. Evidently the physical processes responsible for formation and destruction of PAHs in normal star-forming galaxies tend to maintain PAH abundances near 3% provided that the metallicity 18) it appears that there is a threshold metallicity for PAH formation: 6.6. Dependence of Global Dust-to-gas Ratio on Metallicity 6.6.1. Theoretical Expectations The abundance of dust in the ISM is the result of competition between processes that form dust (dust formation in stellar outflows and dust growth in the ISM) and processes that return material to the gas phase (e.g., sputtering in hot gas and vaporization in high-speed grain-grain collisions). In the Milky Way and other star-forming galaxies with near-solar metallicity, accretion of atoms onto grains is rapid in the cool, dense phases of the ISM, and the balance between grain growth and grain destruction maintains a large fraction of the refractory elements in grains. Most of the dust in the Milky Way must have been grown in the ISM; there is simply no other way to understand the observed extreme depletions of elements like Si, Al, Ca, Ti, and Fe in the diffuse ISM (see, e.g., Draine 1990;Weingartner & Draine 1999;Draine 2009).
The black dashed line in Figure 8(a) shows the expected dependence of M d /M H on O/H if all galaxies had heavy element abundances proportional to solar abundances and the same depletion pattern as measured in the well-studied cloud toward the nearby star ζOph; in this cloud, the refractory elements (e.g., Mg, Si, Fe) are almost completely incorporated into grains, and we infer a dust/H mass ratio of 0.0099 (see Table 23.1 of Draine 2011). For this scenario, we then expect Asplund et al. 2009, corrected for diffusion). However, in the overall ISM, M d /M H will fall below this limiting value, because of dust destruction processes. A simple toy model can illustrate the competition between formation and destruction processes (similar models have been discussed by, e.g., Edmunds 2001;Mattsson et al. 2012;Asano et al. 2013).
Let Z m be the fraction of the ISM mass in "refractory" elements ( »Ź ). Clearly < Z Z d m , since some of the refractory elements are in the gas phase.
Destruction and grain growth in the ISM both contribute to the rate of change of Z d . We also include a term representing injection of solid grains into the ISM from stellar sources (AGB stars, supernovae, and so on). The rate of change of Z d is given by The first term t -Z d d is the rate of dust destruction: τ d is the lifetime of solid material in the ISM against destructive processes that return material to the gas phase. The destruction rate td 1 is a mass-weighted average over the dust in the (see discussion in, e.g., Draine 2009). Realistic estimation of τ d requires a detailed, dynamic multiphase model of the ISM (e.g., Zhukovska et al. 2016). The appropriate value of τ d will obviously vary with galactocentric radius within a galaxy and from galaxy to galaxy. The term in Equation (20) representing grain growth is proportional to ) because it depends on grain surface area (µZ d , for a fixed distribution of grain sizes) and on the gas-phase abundance of condensible elements is the probability per unit time that a refractory atom in the gas phase will collide with and stick to a grain.
The last term, t  Z m , represents injection of dust into the ISM from stellar sources, such as cool AGB stars, planetary nebulae, and core-collapse supernovae. This term will obviously depend on the stellar populations. Here, for illustration, we take the injection rate to be proportional to the metallicity Z m . For the galaxies of interest here, this injection term is small compared to the other terms in Equation (20), and the precise form adopted in Equation (20) is not critical.
If the shortest of the timescales t t , a d { } is short compared to the ∼10 yr 9 timescale for galactic chemical evolution, and , we can neglect the time dependence of the metallicity Z m . The toy model will approach a quasi-steadystate solution with This solution for Z d depends only on Z m and on ratios of timescales, t t a d and t t  : dust formation in stellar outflows is secondary to dust growth in the ISM (i.e., only a small fraction of interstellar dust is "stardust").
For large values of Z m , all models approach the upper limit = Z Z 1 d m (long-dashed line in Figure 8(b)). Models of interest have t t < a d , so that for near-solar abundances, accretion is faster than destruction, and a solarmetallicity ISM can maintain a large fraction of the refractory elements in dust (i.e., ). However, for sufficiently low O/H, accretion rates become slow, resulting in low values of Z d /Z m .

Observations
Using dust mass estimates based on modeling the infrared emission, radial variations in DGRs were found for galaxies in the SINGS sample (Muñoz-Mateos et al. 2009) and for M101 (Vílchez et al. 2019 the dust-to-metals ratio appeared to decline with decreasing metallicity. Chiang et al. (2018) found variations in the dust-to-metals ratio in M101, which they related to variations in both metallicity and H 2 fraction. De Cia et al. (2016) found similar behavior in a sample that included 55 damped Lyα systems (DLAs), where dust abundances were inferred from depletions of Si and metallicities from [Zn/Fe]. It appears that as metallicity decreases below a certain threshold (e.g., ), an increasing fraction of refractory elements (Mg, Si, Fe, ...) remains in the gas phase.
DGRs for the KF57 sample (see Table 1) are plotted against O/ H in Figure 8(b), with dust masses estimated from our model, gas masses taken from Table 10, and the PP04N2 estimate for O/H. Fourteen galaxies have detections of both dust and H I, but were either not observed or not detected in CO, resulting in DGR upper limits. An additional seven galaxies were detected in H I but not in dust, resulting in DGR upper limits. Figure 8(b) shows a clear dependence of dust/gas ratio on metallicity. With some exceptions, the observed dust/H mass ratios for the KF57 sample are in broad agreement with the toy model (Equation (19) (green curve in Figure 8 has a measured dust/H mass ratio several times larger than the "upper limit" Z Z 0.0099( )  (although NGC 1482 is missing CO measurements). It is notable that the ISM appears to have been subject to unusual activity. NGC 1482 shows evidence of a galactic-scale "superwind": the X-ray morphology shows a striking "hourglass" shape emerging from the plane of the disk (Strickland et al. 2004;Vagshette et al. 2012). Interestingly, this galaxy is completely missing H I in its central region, with atomic gas only found in two blobs ∼2 kpc distant from its center, roughly at the confines of the X-ray emission (Hota & Saikia 2005 ( ) ( )  , the M d /M H ratio would be normal for its metallicity. 2. NGC 4594 (M104 "Sombrero," type SAb) also has nearsolar O/H, but a dust/H mass ratio several times larger than the expected upper limit Z Z 0.0099( )  . NGC 4594 has diffuse X-ray emission, suggesting the presence of a galactic-scale outflow (Li et al. 2011). Li et al. (2011 estimate the hot gas to have a temperature »T 6 10 K 6 and total mass »Ḿ M 2.9 10 hot 8  . Adding this to the Bajaja et al. (1984) value for H I and the H 2 mass estimated with a standard X CO factor, we find =Ḿ M 6.0 10 , about a factor of 2.5 above the ratio expected for metallicity » Z Z 1  . The gas in the hot phase, with a density » n 0.1 cm H 3 , has a cooling time t »5 10 yr 7 (Li et al. 2011). Some of the hot gas may have cooled down to ∼10 K 4 , perhaps making an additional contribution to the total gas mass present in NGC 4594. We suggest that NGC 4594 may contain a substantial mass of diffuse H II at ∼10 K 4 that has not yet been detected. Gravity, radiation pressure, and inertia can all lead to velocity differences between gas and dust, allowing the two to separate. However, because dust is generally well coupled to the gas by both gas drag and the Lorentz force on charged grains, gas-grain "slip" velocities are generally small (e.g., Weingartner & Draine 2001b), and scenarios where gas is removed but dust is left behind are not viable unless the gas flows are slow enough that the small  (Li et al. 2011). gas-grain "slip" velocities suffice to prevent the dust grains from leaving the galaxy. Even if gas is stripped or lost in an outflow, we expect the metallicity in the remaining gas (and therefore the upper bound Equation (19) on the dust/mass ratio) to be unaffected. If NGC 1482 and NGC 4594 truly have high dust/gas ratios, then this would appear to require a mechanism for concentrating the dust in part of the gas and removing the dust-poor gas via an outflow or stripping. Alternatively, perhaps the dust/gas ratio is actually normal, but the dust mass has been overestimated because the dust material for some reason has a far-infrared/submillimeter opacity that is significantly larger than found in normal starforming galaxies. The elevated dust/gas mass ratios in NGC 1482 and NGC 4594 require further study. 3. NGC 2841 (Type SAb) and NGC 5055 (M63, Type SAbc): These two galaxies have much lower dust/gas ratios than would be expected given their high estimated metallicities . The photometry for these galaxies is reliable, and the models reproduce the SED out to 500μm. However, it seems likely that the metallicities given in Table 2 are overestimated. Moustakas et al. (2010) with KK04 found that NGC 2841 and NGC 5055 have + 12 log O H 10 ( ) significantly greater than 9.0, ∼0.2 dex higher than the values found for the central regions in the same galaxies by Pilyugin et al. (2014). Pilyugin et al. (2014) also deduce strong metallicity gradients in these two galaxies, implying that at R 0.5 25 , the characteristic + 12 log O H 10 ( )∼8.6. Such metallicities, better representing the average over the galactic disk, would be consistent with the observed dust/gas ratios for these two galaxies.
In Figure 8

Resolved Trends of DL07 Parameters
Using data at M160 resolution, the synergy of Herschel and Spitzer for the KINGFISH sample enables an assessment of dust properties on kiloparsec scales in nearby galaxies (the FWHM of the M160 PSF, 38 8, corresponds to 1.86 kpc at the median KINGFISH sample distance of 9.9 Mpc). The number of M160 ´ 18 18 pixels in each galaxy ranges from 20 for the smallest galaxies (M81 dwB, NGC 584) to >4000 pixels for the largest ones (NGC 5457=M101 and IC 342); the resolved sample as a whole, including the nine "extra" galaxies, comprises >32,000 pixels with well-defined dust parameters and photometry. Figure 9 shows how the starlight intensity parameters U min,DL07 and U DL07 are distributed over the ∼32,000 galaxy mask pixels (i.e., S > S Ld Ld,min ) where we are able to estimate the dust and starlight parameters. Half of the pixels have < U 1, and half of the pixels have < U 1 min . The Ū distribution for the KF62 sample ( Figure 9) is similar to that for Local Group galaxies (Utomo et al. 2019).
The distributions of dust luminosity densities S Ld and dust mass densities S Md are displayed in Figure 10.  , corresponding to a starlight heating rate parameter » U 1. The dust light and mass surface densities that most contribute to the total dust budget are more clearly seen in Figure 11, where we show the cumulative distributions of dust luminosity L d and dust mass M d plotted against S Ld and S Md , respectively.
The vertical dashed lines in Figure 11 show the surfacedensity thresholds that provide 50% of the total: S = . Regions with dust light and mass surface densities greater than these values comprise only a small fraction of the total: from Figure 10 we see that 50% of the dust light comes from only ∼3% of the (brightest) pixels, and 50% of the total dust mass from ∼22% of the (densest) pixels.
In what follows we have applied a limit in dust surface brightness S´- L 2 10 kpc ; Ld 6 2  thus the low-S/N, faint outer regions of the sample galaxies (where estimates of parameters such as S Md and q PAH may become unreliable) will not be considered. As seen above, such regions contribute very little to either the light budget or the mass budget of the dust Figure 9. Distributions of U min and U (Ubar) for all galaxies. U min and U are the renormalized values (see Equations (13)-(16)). Both distributions are fairly broad; for a given pixel, U -U min may be smaller than or comparable to U min given that f PDR is usually modest. Thus, it is not surprising that the two peaks are similar. 41 Using PT metallicities, Rémy-Ruyer et al. (2014)  over the sample as a whole. Applying such a cut ensures that the plotted DL07 parameters (and the photometric quantities) will be as accurate as possible, given the constraints of the data; the total number of ´ 18 18 pixels in the sample is reduced to ∼25,500.
We now investigate the IR observational signatures associated with dust heating (U min ). Figure 12 shows U min for all galaxies plotted versus MIPS and SPIRE flux density ratios, f f 70 160 , f f 70 250 , and f f 160 500 . Because of the unexplained discrepancies between MIPS and PACS photometry (see Figure 2), we have elected to use only MIPS photometry for f 70 and f 160 . The left panel shows that the flux ratio f f 70 160 is not a very good predictor of U min . This is because when  U 1 min , the 70 μm emission has an appreciable contribution from (1) single-photon heating of small grains and (2) dust in regions with high starlight intensities (assuming g > 0, which is almost always the case). The flux ratio f f 160 250 , shown in the middle panel, ameliorates the potential domination of the emission by small-grain stochastic heating, but the wavelength ratio of the two fluxes is insufficient to reliably sample U ; min a small range in flux ratio corresponds to as much as an order of magnitude change in U min . However, the right panel shows that the f f 160 500 flux ratio correlates quite well with U min because the emission at both 160 and m 500 m is dominated by the larger grains heated by starlight intensities near U min . Because 160 μm is not in the Rayleigh-Jeans limit for the grain temperatures in these galaxies, the f f 160 500 ratio is sensitive to large-grain temperature, and hence to starlight heating rate. The best-fit correlation, obtained with median clipping and a "robust" regression algorithm, effective for minimizing the effects of outliers (R Core Team 2014), is given by This relation predicts U min to within 0.21 dex (rms) over a range of U min of more than two orders of magnitude. Because the emission at these wavelengths is dominated completely by large grains, this long-wavelength ratio predicts very well the minimum starlight heating intensity. The PAH abundance parameter q PAH varies from galaxy to galaxy, as discussed in Section 6.5, where it is apparent that there is a correlation between q PAH and the gas-phase metallicity O/H. Note that q PAH also exhibits significant variations within individual galaxies, as can be seen from the map of q PAH in M101 (see Figure 5), as well as for other wellresolved galaxies (see ). If q PAH is sensitive to metallicity, then we may expect radial variations within galaxies, with q PAH generally declining with radius. However, our q PAH maps also exhibit substantial azimuthal variations, suggesting that the PAH abundance responds to changes in environmental conditions beyond metallicity alone.
In Figure 13, we explore-using three different proxies for the starlight intensity-whether q PAH is affected by the intensity of the radiation field. The left panel in Figure 13 indicates that q PAH seems to be relatively independent of variations in the f f 70 160 flux ratio. The f f 70 160 flux ratio is apparently not uniquely tracing the temperature of the larger grains; as seen in Figure 12 and discussed below, this ratio begins to reflect U min , and thus large-grain temperature, only above a certain U min threshold ( ) , but the right panel shows a stronger trend where q PAH tends to  ) arises because single-photon heating generally dominates at 24 mm; big grains only get hot enough to radiate at 24 μm when the radiation field is extremely intense. Instead, at 70 μm, singlephoton heating makes a significant contribution only for is a better indicator of warm, large grains than L L 24 d ( ) , and it is these warm, large grains that are the signature of high-intensity radiation fields that could be associated with PAH destruction.
As seen in the right panel of Figure 13, the PAH fraction appears to vary with L L 70 Such a trend may reflect a tendency for PAH destruction to occur in star-forming regions, where O stars supply highenergy photons that photodestroy PAHs, and a significant fraction of the dust is exposed to starlight intensities high enough to elevate the L L 70 dust ( ) ratio. Many studies have previously noted suppression of PAH emission in H II regions is not a good indicator of U min , because f 70 is sensitive to both single-photon heating and the emission from dust exposed to starlight intensities > U U min . The middle panel with f f 160 250 avoids using f 70 , but the wavelength range is insufficient to adequately sample U min and a luminosity-weighted dust temperature. Instead, the right panel shows the tight correlation between U min and f f 160 500 (the dashed line is Equation (22)), illustrating the close relationship between the minimum heating intensity and the coolest dust. (e.g., Giard et al. 1994;Helou et al. 2004;Povich et al. 2007;Relaño et al. 2016). In a detailed study of PAH abundances in the Magellanic Clouds, Chastenet et al. (2019) show that q PAH is reduced in regions close to sources of H-ionizing radiation.
High values of L L 70 d ( ) also occur in low-metallicity galaxies, because of radiative transfer effects (hotter stars, less dust attenuation), and is consistent with the tendency for lower q PAH in a metal-poor ISM.

Resolved Dust Light-to-mass Ratios
The dust light-to-mass ratio in galaxies and within galaxies, L d /M d , should reflect the peak and spread of luminosityweighted dust temperatures. In the DL07 model, µ U L d /M d , so U also probes dust temperatures. We would thus expect L d /M d to depend on photometric flux ratios, as long as the two wavelengths in the flux ratios are sampling a sufficiently broad spectral range to be sensitive to large-grain temperature variations. Figure 14 illustrates the correlations in the resolved pixels of all galaxies between S Ld /S Md (noted as L d /M d in the ordinate axis label) and, as in Figure 12, three flux density ratios, f f 70 160 , f f 160 250 , and f f 160 500 . The longer the wavelength ratio (in this case 160 μm/500 μm), the better that S Ld /S Md can be predicted from observations. The right panel (black long-dashed line) of Figure 14 shows the correlation with f f 160 500 given by This fit with f f 160 500 has an rms deviation of 0.16 dex over >22,000 degrees of freedom (dof). The trend of L d /M d with f f 160 250 is much less reliable, so we have not shown any regression in the middle panel of Figure 14. Because of the limited wavelength lever arm for the f f 160 250 flux ratio (see also Figure 12), for a given f f 160 250 ratio, L d /M d can vary by a factor of 30 or more; this makes it difficult to accurately determine the dust light-to-mass ratio from f f 160 250 . Nevertheless, if we know the dust luminosity L d and have a measure of a flux around the peak of dust emission (e.g., f 160 ) and one sufficiently far away and in the Rayleigh-Jeans regime (e.g., f 500 ), we can estimate the dust mass M d to within ∼50%. Figure 15 shows the same quantities but separately for three galaxies representative of the extremes probed by the KING-FISH sample: IC 2574, a metal-poor dwarf; NGC 5457 (M101), a face-on grand-design spiral; and NGC 2146, a luminous IR galaxy (LIRG). For a flux density ratio with short +long wavelengths (e.g., f f 70 160 ), the L d /M d ratio within these galaxies can differ by up to an order of magnitude. As has been seen in previous figures, because f f 70 160 is sensitive to both single-photon heating and the possible exposure of a small fraction of the dust to starlight intensities > U U min , the f f 70 160 ratio does not strongly constrain the temperature of the dust grains that dominate the total emission. Instead, the longer wavelength ratio (e.g., f f 160 500 ) is a much better indicator of large-grain temperature, and consequently better correlated with the dust light-to-mass ratio L d /M d . Two regressions are shown in Figures 14 and 15; the (black) long-dashed line, described above (see Equation (25)), is for the entire sample. The (gray) dashed-dotted one is the regression obtained for only IC 2574 and NGC 2146 and is given by The regression for the entire sample is entirely consistent with NGC 5457 (M101), but not for IC 2574 and NGC 2146, which may be considered two "extreme" galaxies. The overall radiation fields U in IC 2574 and NGC 2146 are higher (in the mean, by ∼40% and a factor of 13, respectively) than that of NGC 5457. These more intense heating fields, possibly a signature of starbursts, result in a slightly steeper slope relating L d /M d and f f 160 500 than in more quiescent environments such as the disk of NGC 5457 (and most of the KINGFISH sample).
In the present model, the dust temperatures are determined by the starlight intensity distribution within a pixel, which is , and f f 160 500 (right) for all galaxies. The color coding corresponds to pixel number density, as shown by the rightmost color bar. The best-fit (robust) regression for f 160 /f 500 is shown as a (black) long-dashed line and corresponds to rms residuals of ∼0.18 dex (see Equation (25)). The (gray) dashed-dotted line is the analogous best-fit regression for only IC 2574 and NGC 2146 (with 383 dof; see Figure 15). characterized by three parameters, U min , γ, and α (see Section 5), and we need more than two bands if we wish to determine the distribution of temperatures for the emitting dust well enough to reliably estimate the mass of dust in the pixel. On the other hand, if flux ratios at longer wavelengths are considered (right panel of Figure 15), there is much less variation within and between galaxies. Such behavior was also seen with U min in Figure 12 and suggests that ratios at these longer wavelengths better trace L d /M d because they provide better information about the temperatures of the large grains that dominate the dust luminosity.

Resolved Dust Mass Surface Densities
Given the relative constancy of dust-to-metals ratios for galaxies with metallicities +  12 log O H 8.4 10 ( ) (see Figure 8(b)), the dust luminosity in the Rayleigh-Jeans (R-J) regime of the dust SED has become a popular tracer of ISM mass (e.g., Corbelli et al. 2012;Eales et al. 2012;Scoville et al. 2014Scoville et al. , 2016Scoville et al. , 2017Groves et al. 2015). This is an effective technique both locally and at high redshift because the R-J tail of the dust emission probes optically thin dust and is relatively insensitive to dust temperature. Here we explore whether this is also true for the spatially resolved dust emission in the KINGFISH sample. Figure 16 shows the dust mass surface density S M d (estimated from the renormalized DL07 model) plotted against monochromatic dust luminosity surface density pn S = n n n I 4 L in the MIPS 160 μm, SPIRE 250 μm, and SPIRE 500 μm bands. It can be seen that at 160 μm, a wavelength that generally probes the dust emission peak, there is only a broad correlation with more than an order of magnitude dispersion at low surface brightness. As wavelength increases toward the SPIRE bands, the correlation improves and becomes very good at 500 μm, similar to the trends found for KINGFISH global values by Groves et al. (2015).
The rightmost panel reports the best-fit correlation, obtained with the robust regression algorithm: , and f 160 /f 500 (right) for three galaxies separately: IC 2547, a lowmetallicity dwarf; NGC 5457 (M101), a large grand-design spiral; and NGC 2146, an LIRG. The contours reflect the individual galaxies (IC 2547 blue, NGC 5457 green, and NGC 2146 red) and correspond to pixel number densities. In the right panel, the (black) long-dashed line corresponds to the best-fit regression reported in Figure 14 for the sample as a whole (see Equation (25)), and the (gray) dashed-dotted line corresponds to the analogous best-fit regression for only IC 2574 and NGC 2146 (383 dof). , and SPIRE 500μm (right) bands, for all galaxies. The color coding corresponds to pixel number density, as shown by the rightmost color bar. The dashed line in the right panel (SPIRE500) shows the best-fit regression relating dust mass surface density S Mdust to 500μm luminosity surface density n m S n 500 m L ( ) (Equation (27)). See Section 6.7.2 for more details.
This fit gives an rms scatter σ=0.07 dex on log 10 (M d ) (with ∼25,400 dof), implying that dust mass surface densities can be inferred from 500 μm luminosity surface densities to within ∼20%. The slope is significantly sublinear, over almost three decades of 500 μm luminosity surface densities, reflecting the tendency for dust to be somewhat warmer in pixels where S Md is high, presumably because these pixels are more likely to harbor star-forming regions. Groves et al. (2015) obtained a similar result globally for inferring gas mass from L 500 for all KINGFISH galaxies, including dwarfs (stellar mass 10 9 M  ); however, once Groves et al. (2015) considered only the more massive galaxies, the slope steepened and became approximately linear.
The rms deviation of only 0.07 dex from Equation (27)  . To estimate ISM mass from Equation (27), the dust mass from Equation (27) needs to be combined with a gas-to-dust ratio, as discussed in Section 6.6. However, this ratio depends on metallicity (see Figure 8); thus, oxygen abundance needs to be incorporated to estimate gas mass for metal-poor galaxies. In any case, Figure 16 shows that the slope between dust mass and luminosity is steeper, closer to unity, at lower surface brightnesses, roughly independent of wavelength. However, global integrated values of quantities such as long-wavelength IR luminosity are luminosity weighted, thus sampling preferentially higher surface brightnesses. Thus, our new result for resolved regions in KINGFISH galaxies is inconsistent with a strictly linear trend of dust mass with long-wavelength IR luminosity. Indeed, as noted above, a nonlinear behavior would be expected since the dust in high S Md pixels is, on average, somewhat warmer.

Summary
Dust modeling results for 70 galaxies (61 KINGFISH galaxies, plus nine additional galaxies present in the observed fields) are presented here. Dust is detected reliably in 62 galaxies, and upper limits are reported for the remaining eight. Tables 5 and 6 report the global galaxy photometry, and the best-fit dust parameter estimates are given in Table 9. Dust parameter maps are displayed in . The DL07 dust model successfully reproduces the dust SEDs over the wide variety of environments present in the KINGFISH sample.
Long-wavelength imaging can be omitted in order to increase the angular resolution of the modeling, but results become unreliable if the long-wavelength coverage is insufficient. For maximum reliability, we recommend using all cameras available, including MIPS160, SPIRE250, SPIRE350, and SPIRE500. If better angular resolution is critical, the lowest-resolution cameras (SPIRE500 and MIP160) can be left out, but estimates of dust mass become unreliable unless at least SPIRE250 is included. If SPIRE350, SPIRE500, and MIPS160 are not included, the DL07 model dust masses can be low by as much as a factor of 0.8 or high by as much as a factor of 2 (see Figures 20); the median factor is 1.25. The q PAH and f PDR estimates are fairly insensitive to the camera combination used, so they can be obtained reliably without l m > 250 m photometry, provided that the S/N is adequate.
Resolved (multipixel) modeling and global (single-pixel) modeling generate similar estimates of M d , q PAH , and f PDR when all the Spitzer and Herschel cameras are employed. The single-pixel modeling tends to slightly underestimate the total dust mass M d by ∼13% (see Figure 21).
Our analysis shows that q PAH , the fraction of the dust mass contributed by PAHs, correlates much better with the PP04N2 estimate for O/H than for the PT estimate, strongly suggesting that PP04N2 is a better strong-line abundance estimator than the PT estimator. We find that q PAH appears to increase monotonically with increasing metallicity, with q PAH varying linearly with log O H ( ) for + > 12 log O H 7.94 10 PP04N2 ( ) (see Figure 7(b) and Equation (18)).
For most star-forming galaxies with metallicity  Z Z  , the dust/gas ratio is close to the limiting value where nearly all of the refractory elements are locked up in grains. However, at lower metallicity, the dust/gas ratio is often well below this limiting value, consistent with what is expected from a simple toy model with accretion rate t µ -Z a d 1 (see Figure 8(b)). The resolved regions in the KINGFISH galaxy sample show several trends with U min , q PAH , and mass-to-light ratios for dust emission. The characteristic starlight intensity U min can be estimated from long-wavelength flux ratios (e.g., f f 160 500 ) to within a factor of two over more than two orders of magnitude in U min (see Equation (22)). From the same flux ratio, and with a measurement of dust luminosity, dust mass can be estimated to within ∼50% (see Equation (25)). Despite a variation of 3 orders of magnitude in IR surface brightness, for the adopted physical dust model it is possible to estimate dust mass from IR luminosity at 500 μm to within ∼0.07 dex), affording an accuracy of ∼20% (see Equation (27)). There are of course systematic errors coming from the choice of dust model, but these are difficult to estimate. Estimating gas mass for metalpoor galaxies requires incorporating metallicity, because of the metallicity dependence of DGRs. Our formulations for inferring starlight heating intensity and dust mass from flux ratios and integrated IR or monochromatic luminosities have been calibrated over 22,000 independent regions in 62 galaxies, spanning metal-poor dwarf irregulars to grand-design spiral disks and actively star-forming LIRGs. These calibrated prescriptions are designed with the aim of facilitating comparison with high-redshift galaxies, where frequently restframe f 160 and at least one longer wavelength flux are available. Figure 17. Hol2: Model results at M160 PSF (rows 1 and 2) and at S250 PSF (rows 3 and 4). Dust luminosity per area S Ld (column 1, rows 1 and 3) is shown for the entire field, with the adopted galaxy mask boundary in white. Dust mass per area S Md (column 2, rows 1 and 3) is after renormalization (see text). U min,DL07 , q PAH , and f PDR are shown in rows 2 and 4. The global SED (column 3, rows 1 and 3) is shown for single-pixel modeling, with contributions from dust heated by U min (green) and dust heated by > U U min (red) and starlight (cyan); values of U min and M d in the figure label are for the DL07 model before renormalization. Herschel (blue rectangles) and Spitzer (red rectangles) photometry is shown; the vertical extent is ±1σ. Diamonds show the band-convolved flux for the model.

Appendix A Resolved Dust Parameter Maps for KINGFISH Galaxies
As described in the text, each galaxy where we have a positive dust detection has two figures: the first (a) shows the model done at MIPS160 resolution, using data from all cameras (IRAC, MIPS, PACS, and SPIRE). This is our "gold standard" modeling. The second (b) shows a model at SPIRE250 resolution, using IRAC, MIPS24, PACS, and SPIRE250 cameras (i.e., omitting MIPS70, MIPS160, SPIRE350, and SPIRE500). This latter modeling, while able to resolve smaller scale structures in the galaxies, is overall less reliable.
Each figure in the Figure 17 figure set has 12 panels. For each of the resolutions, the top row is a map of dust luminosity surface density S L d (left), dust surface density S M d (center), and the model SED (right). The lower row shows the starlight intensity parameter U min,DL07 (left), the PAH abundance parameter q PAH (center), and the PDR fraction f PDR (left). The dust luminosity surface density S L d is shown for the full field, with the white contour showing the minimum surface brightness S Ld,min below which we do not attempt to model the emission. Maps of derived quantities (S M d , U min , q PAH , and f PDR ) are limited to the "galaxy mask" region with S > S L Ld,min d . In the SED plot, the observed photometry is represented by rectangular boxes (Spitzer IRAC and MIPS in red; Herschel PACS and SPIRE in blue) showing ±1σ uncertainties. The black line is a single-pixel DL07 model that seeks to reproduce the observed SED, with different components shown. The values of U min and M d in the label are for the DL07 model before renormalization. The cyan line is the stellar contribution, the dark red line is the emission from dust heated by the power-law U distribution, and the dark green line is emission from dust heated by = U U min .

Appendix B Cases with 3σ Upper Limits for Dust Mass
Eight galaxies in the KINGFISH sample, five dwarfs (DDO 053, DDO 154, DDO 165, Hol1, and M81dwB) and three ellipticals (NGC 0584, NGC 0855, and NGC 1404), yield upper limits on the dust mass from our modeling. In these cases, the signal from the galaxy in the far-IR is of comparable magnitude to the contamination from background galaxies, leading to uncertain dust mass measurements.
For the three elliptical galaxies, we employ a S Ld -based mask that roughly coincides with the optical galaxy. To obtain an upper limit on the dust mass in the case of these nondetections, we randomly shift the mask around in the M160-resolution dust mass image, avoiding overlap with the original mask, and remeasure the total dust mass. We construct a distribution of these "background" dust mass measurements to obtain a mean and standard deviation. In all galaxies mentioned above, the mean of the "background" dust mass values is positive, as would be expected due to the real signal at far-IR wavelengths from unidentified background sources ("confusion noise"). The measured dust mass at the expected galaxy location is within ∼1σ of the mean. In the text and Table 11, we provide the 3σ upper limit on the dust mass generated with this procedure.
The definition of the galaxy mask itself is also potentially affected by confusion noise. In the case of the dwarf galaxies, we use the H I observations from THINGS and LittleTHINGS to create an alternative galaxy mask, based on a cut at an H I column density of -10 cm 20 2 from the H I image convolved to M160 resolution. The H I-based galaxy mask is typically somewhat larger than that defined by the dust luminosity surface density cut. We apply the same procedure described above to obtain the background mean and standard deviation. Table 11 and Figure 18 provide the results of this procedure. Table 11 lists the 3σ upper limits for each galaxy using the H Ibased masks for the dwarfs and the dust luminosity surface density masks (as described in the text) for the ellipticals. Histograms of the dust masses from the randomly shifted masks are shown in Figures 18.1-18.8. We note that the dust mass limits from this procedure are expected to be very conservative. Higher S/N could be obtained by a careful treatment of the integrated photometry for each galaxy, taking the confusion noise into account.  Figure 19 shows the comparison of the dust parameter estimates obtained from models using the S250 PSF (18.2″ FWHM) with parameters estimated from (our "gold standard") modeling using all cameras (IRAC, MIPS, PACS, SPIRE) and the M160 PSF (39″ FWHM).
The results for L d and M d at S250 resolution appear to be quite robust: the median change in L d is only 5%, which may be due in part to calibration differences between MIPS160 and PACS160, with MIPS160 data being used only in the M160 PSF modeling. Note that M d shows more variation, with a median change of 25%; this is likely because loss of SPIRE350 and SPIRE500 may allow modeling at the S250 PSF to include a bit more cool dust than is actually present. However, it is gratifying that the median change is only 25%, indicating that the DL07 model is relatively good at "predicting" l m > 300 m emission using data shortward of m 300 m. However, in some cases, the dust mass is overestimated by as much as a factor of 2 (see Figure 20), and we therefore recommend using M160resolution modeling rather than the riskier S250 PSF.  This is the riskiest PSF we are willing to consider.

D.2. Dust Mass Estimates at Different Resolutions
For each resolution, Figure 20 shows a histogram of the galactic total dust mass estimates divided by the gold standard estimate. We observe that dust mass discrepancies can be large, with the errors and bias increasing as fewer cameras are used, and long-wavelength data are lost. The S500 case (coverage out to m 500 m, a PSF that is not much smaller than the M160 PSF, but no MIPS160 photometry) gives dust mass estimates that are close to our gold standard estimate, with a median ratio of 1.21. However, there are a few outliers where M d appears to be overestimated by as much as a factor of 2. These are all galaxies with very weak dust emission and low-S/N data, where loss of the data from one camera (MIPS160) causes a significant change in the apparent SED.
The systematic bias in M d and the scatter both increase as we move to smaller PSFs (S350, S250, P160). At P160 resolution, fully 25% of the cases have M d under estimated by a factor of 2 or more.
On balance, it appears that modeling at S250 resolution is reasonable, although slightly risky: there is a significant chance that the dust mass may be overestimated or underestimated by a factor of 1.5 or more. S350 resolution is safer, and S500 even better. Figure 19. Comparison of modeling with the S250 PSF to results obtained with the M160 PSF (the "gold standard"). With S250 modeling, the dust mass tends to be slightly overestimated; the median overestimation is ∼25%.

D.3. Modeling Using Only Global Photometry (Single Pixel)
The KINGFISH galaxies are close ( < D 30 Mpc) and large enough that they can be resolved using the Herschel Space Telescope. When studying galaxies at larger distances, only their global photometry may be available. Here we compare our dust mass estimates using resolved imaging and multipixel modeling with "single-pixel" modeling that makes use of only the global SED. We recall that the dust modeling is not a linear process, and differences in parameter estimates are to be expected. Figure 21 shows the ratio of the dust model parameter estimates from fitting global photometry (a "single pixel" model) versus our "gold standard" multipixel modeling at M160 resolution, where each pixel is modeled separately. In both cases we use all cameras: IRAC, MIPS, PACS, and SPIRE. We observe that L d is very reproducible, with the single-pixel luminosity estimate differing from the multipixel result by only a few percent.
The dust mass estimate is probably most important, and is found to be moderately robust: for 75% of the cases, the singlepixel modeling obtains a mass estimate within 25% of the resolved multipixel analysis. Thus, dust mass estimates for unresolved distant galaxies should be reliable, assuming only that the photometry covers a suitable range of wavelengths (rest-frame wavelengths l m   50 300 m), with an adequate S/N.
The á ñ f PDR and á ñ q PAH estimates are both quite robust, with the single-pixel results agreeing with the multipixel analysis to within ∼5% in most cases. Figure 20. Comparison of estimates for dust mass M d obtained from multipixel modeling with different PSFs and camera combinations. The "gold standard" refers to modeling with the M160 PSF and all cameras. Here we compare M d obtained with the S500, S350, S250, and P160 PSFs (see Table 4). Because of the limited wavelength coverage and lower S/N, for the P160 PSF the dust mass can be in error by up to a large factor: 10/62 cases underestimate the dust mass by more than a factor of 2, and 2/62 cases by more than a factor of 5.