# The Resolved Behavior of Dust Mass, Polycyclic Aromatic Hydrocarbon Fraction, and Radiation Field in ~800 Nearby Galaxies

Jérémy Chastenet, Karin Sandstrom, Adam K Leroy, Caroline Bot, I-da Chiang, Ryan Chown, Karl D Gordon, Eric W Koch, Hélène Roussel, Jessica Sutter, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jérémy Chastenet, Karin Sandstrom, Adam K Leroy, Caroline Bot, I-da Chiang, et al. The Resolved Behavior of Dust Mass, Polycyclic Aromatic Hydrocarbon Fraction, and Radiation Field in ~800 Nearby Galaxies. The Astrophysical Journal Supplement Series, 2024, 276 (1), pp.2. 10.3847/1538-4365/ad8a5c. insu-04878751

HAL Id: insu-04878751

https://insu.hal.science/insu-04878751v1

Submitted on 10 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The Resolved Behavior of Dust Mass, Polycyclic Aromatic Hydrocarbon Fraction, and Radiation Field in ~800 Nearby Galaxies

Jérémy Chastenet¹, Karin Sandstrom², Adam K. Leroy³, Caroline Bot⁴, I-Da Chiang (江宜達)⁵, Ryan Chown⁶, Karl D. Gordon¹,⁷, Eric W. Koch⁸, Hélène Roussel⁹, Jessica Sutter¹⁰, and Thomas G. Williams¹¹

¹Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, B-9000 Gent, Belgium; jeremy.chastenet@ugent.be

²Department of Astronomy & Astrophysics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA

³Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 West 18th Avenue, Columbus, OH 43210, USA

⁴Observatoire Astronomique de Strasbourg, Université de Strasbourg, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France

⁵Institute of Astronomy and Astrophysics, Academia Sinica, No. 1, Sec. 4, Roosevelt Road, Taipei 106216, Taiwan

⁶Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, OH 43210, USA

⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

⁸Center for Astrophysics—Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

⁹Sorbonne Université, Institut d'Astrophysique de Paris, CNRS (UMR7095), 75014 Paris, France

¹⁰Whitman College, 345 Boyer Avenue, Walla Walla, WA 99362, USA

¹¹Subdepartment of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
Received 2024 July 12; revised 2024 September 16; accepted 2024 September 20; published 2024 December 11

## Abstract

We present resolved 3.6–250 µm dust spectral energy distribution (SED) fitting for ~800 nearby galaxies. We measure the distribution of radiation field intensities heating the dust, the dust mass surface density ($\Sigma_d$), and the fraction of dust in the form of polycyclic aromatic hydrocarbons (PAHs; $q_{PAH}$). We find that the average interstellar radiation field ($\bar{U}$) is correlated both with stellar mass surface density ($\Sigma_*$) and star formation rate surface density ($\Sigma_{SFR}$), while more intense radiation fields are only correlated with $\Sigma_{SFR}$. We show that $q_{PAH}$ is a steeply decreasing function of $\Sigma_{SFR}$, likely reflecting PAH destruction in H II regions. Galaxy-integrated $q_{PAH}$ is strongly, negatively correlated with specific star formation rate (sSFR) and offset from the star-forming “main sequence” (ΔMS), suggesting that both metallicity and star formation intensity play a role in setting the global $q_{PAH}$. We also find a nearly constant $M_d/M_*$ ratio for galaxies on the main sequence, with a lower ratio for more quiescent galaxies, likely due to their lower gas fractions. From these results, we construct prescriptions to estimate the radiation field distribution in both integrated and resolved galaxies. We test these prescriptions by comparing our predicted $\bar{U}$ to results of SED fitting for stacked “main-sequence” galaxies at $0 < z < 4$ from M. Béthermin et al. and find sSFR is an accurate predictor of $\bar{U}$ even at these high redshifts. Finally, we describe the public delivery of matched-resolution Wide-field Infrared Survey Explorer and Herschel maps along with the resolved dust SED-fitting results through the Infrared Science Archive.

*Unified Astronomy Thesaurus concepts: Interstellar dust (836); Polycyclic aromatic hydrocarbons (1280); Infrared photometry (792)*

## 1. Introduction

Interstellar dust is fundamentally important to the physics of the interstellar medium (ISM). By absorbing photons over a wide range of wavelengths and converting their energy to grain heating (which is then reradiated primarily in the infrared) or gas heating (via the ejection of photoelectrons), dust couples starlight to the ISM. The reradiated starlight also serves as a tracer of star formation observable across a wide redshift range (R. C. Kennicutt & N. J. Evans 2012; L. J. Tacconi et al. 2020). The amount, composition, and size of grains govern their absorption and scattering efficiencies, infrared emissivity, photoelectric yields, and other key properties (e.g., B. T. Draine 2003; F. Galliano et al. 2018; B. S. Hensley & B. T. Draine 2023; N. Ysard et al. 2024). It is therefore of great importance for a wide range of astrophysical topics to understand what sets the abundance and properties of dust and the radiation field illuminating it.

Infrared (IR) emission is an important tool to study dust because it represents the most accessible direct tracer of dust, and it can be observed across a wide range of galaxy types and redshifts. Extinction and elemental depletions provide crucial complementary constraints on grain-size distributions and composition (for a recent compilation of such constraints, see B. S. Hensley & B. T. Draine 2021, and references therein). However, in galaxies outside the Local Group, the necessary measurements to study depletion and extinction can only be done in small samples of sight lines with sufficiently bright background point sources (e.g., bright stars in relatively nearby galaxies, or background quasars as sample dust in their foreground; E. B. Jenkins 2009; C. Péroux & J. C. Howk 2020; J. Roman-Duval et al. 2022). Measurements of dust attenuation can also provide insights into dust content and properties, but complexities related to geometry and scattering make such studies challenging (for a review, see S. Salim & D. Narayanan 2020). By contrast, dust emission is detectable from the near-IR through millimeter and has been observed extensively by ground, stratospheric, and space-based facilities.

In particular, the Herschel Space Observatory (G. L. Pilbratt et al. 2010), over the course of its mission, observed far-IR dust emission in a large sample of resolved, nearby galaxies, at 70,

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
100, and 160 µm using the Photodetector Array Camera and Spectrometer (PACS; A. Poglitsch et al. 2010) and at 250, 350, and 500 µm with the Spectral and Photometric Imaging REceiver (SPIRE; M. J. Griffin et al. 2010). Herschel observations will remain the best-quality (in sensitivity and angular resolution) far-IR data set until the next-generation far-IR facility comes into service. The set of nearby galaxies observed with Herschel PACS and SPIRE are therefore a critical resource for understanding dust properties, and are the focus of this study.

The basic features of the IR spectral energy distribution (SED) of dust emission are set by (1) the dust mass surface density ($\Sigma_d$), (2) the dust composition (including the fraction of dust in the form of polycyclic aromatic hydrocarbons, PAHs; L. J. Allamandola et al. 1989; see also reviews from A. G. G. M. Tielens 2008; A. Li 2020), and (3) the intensity (and potentially the spectrum) of the radiation field illuminating the dust.<fnref n="12" /> The intensity of the radiation field sets the steady-state temperature of the grains in equilibrium with the radiation field. This temperature can be inferred from the modified blackbody-like emission that the equilibrium grains produce (peaking typically between 100 and 200 µm). Given the equilibrium temperature and knowledge of the dust grain optical properties, the dust mass surface density, $\Sigma_d$, can be inferred from the intensity of the equilibrium IR emission. Small grains are not in equilibrium with the radiation field, but rather are stochastically heated to high temperatures by the absorption of single UV photons (K. Sellgren 1984; B. T. Draine & N. Anderson 1985; B. T. Draine & A. Li 2001). A key component of the small grain population is PAHs, which produce bright vibrational emission bands in the mid-IR after absorbing a UV/optical photon. Given the intensity of the illuminating radiation field constrained by the equilibrium grains, the mid-IR emission allows one to infer the fraction of the dust mass in the form of PAHs (B. T. Draine et al. 2007) and the distribution of radiation field intensities.

Within the beam of typical extragalactic IR observations (hundreds of parsecs to kiloparsec scales), a range of radiation field intensities will inevitably be present (G. Helou 1986; D. A. Dale et al. 2001; D. A. Dale & G. Helou 2002). It is therefore necessary to make some assumptions about the subresolution radiation field distribution in modeling the IR SED. This is particularly important because of the very steep dependence of dust luminosity on temperature, which can lead to biased temperatures (e.g., D. Utomo et al. 2019). One common approach, based on work by D. A. Dale et al. (2001) and further expanded upon by B. T. Draine et al. (2007) involves a power-law distribution of radiation field intensities, with a delta function at a minimum intensity ($U_{min}$) meant to represent the overall diffuse ISM interstellar radiation field. The assumed distribution of radiation field intensities is a critical aspect of interpreting the dust SED. In models with a single radiation field, the width of the equilibrium dust peak and mid-IR emission can be interpreted as indicators of changing dust properties, particularly variations in very small grain abundance or dust composition (see Figure 3 of F. Galliano et al. 2018 for a clear visualization).

A variety of studies have performed IR SED fitting in nearby galaxies, both resolved (J. C. Muñoz-Mateos et al. 2009;

G. Aniano et al. 2012, 2020; K. D. Gordon et al. 2014; L. K. Hunt et al. 2015; V. Casasola et al. 2017, 2022; J. Chastenet et al. 2021; C. J. R. Clark et al. 2021, 2023; Abdurro'uf et al. 2022) and unresolved (A. Boselli et al. 2010; L. Ciesla et al. 2014; A. Rémy-Ruyer et al. 2014; J. I. Davies et al. 2017; C. J. R. Clark et al. 2018; F. Galliano et al. 2021; D. A. Dale et al. 2023), using various dust SED models. These studies have revealed the existence of scaling relationships between dust mass, radiation field, and PAH fraction with galaxy properties and examined the correlations of the parameters with the local galactic environment. Key trends that have been identified include the scaling of dust mass with stellar mass and metallicity (E. da Cunha et al. 2010; L. Cortese et al. 2012; F. Calura et al. 2017; F. Galliano et al. 2021); PAH fraction with stellar mass, metallicity, and radiation field intensity (A. Nersesian et al. 2019; G. Aniano et al. 2020; F. Galliano 2022); and radiation field parameters with star formation rate (SFR) and stellar mass (A. Nersesian et al. 2019). In general, samples where resolved dust SED fitting have been performed are often small, focusing on ~tens of targets (e.g., G. Aniano et al. 2020; V. Casasola et al. 2022). The samples studied with integrated SED fits are typically much larger, up to ~900 galaxies as in the Dustpedia<fnref n="13" /> compilation (J. I. Davies et al. 2017, 2019; S. Bianchi et al. 2018; C. J. R. Clark et al. 2018, 2019; A. Nersesian et al. 2019; V. Casasola et al. 2020; F. Galliano et al. 2021).

In the following work, we extend resolved analysis to ~800 nearby galaxies with observations from Herschel. This study focuses on deriving resolved dust properties in a large sample of nearby galaxies, and the relationship between dust and stellar parameters, such as stellar mass and SFR surface density. Covering a wide range of both dust and stellar properties is essential to fully appreciate the interplay between ISM properties and galaxy evolution. The archive of the Herschel space mission is a goldmine for a resolved study of dust property variations in the $z=0$ Universe.

The paper is organized as follows. In Section 2 we describe the definition of the Herschel archival sample, which includes all nearby galaxies with PACS or SPIRE observations. In Section 3 we describe the reduction of the archival Herschel data, and the additional data processing we perform. We present the details of the fitting procedure in Section 4, and a description of the results in Section 5. In Section 6, we discuss the new insights into the properties of dust and radiation fields across the local galaxy population. Our conclusions are summarized in Section 7, and we describe our data delivery to the NASA Infrared Science Archive (IRSA) in Appendix A (doi:10.26131/IRSA581).

## 2. Galaxy Sample

Our galaxy sample is defined as the overlap between the nearby galaxies in the HyperLeda<fnref n="14" /> database (D. Makarov et al. 2014) and those that are found to be covered by either a PACS or SPIRE photometric observation in the Herschel Science Archive.<fnref n="15" /> We begin with the full sample of galaxies at $cz < 5000 \text{ km s}^{-1}$ in HyperLeda. To create this selection, we query the database for all galaxies (i.e., object type “G”) with measured heliocentric radial velocities <5000 km s⁻¹. This

<footnote n="12">¹² We assume optically thin emission in the mid- to far-IR throughout this study. In the nearby galaxies we study, significant optical depth in the IR is rare. Attenuation can be an issue for ultraluminous infrared galaxies.</footnote>

<footnote n="13">¹³ http://dustpedia.astro.noa.gr</footnote>

<footnote n="14">¹⁴ http://leda.univ-lyon1.fr/</footnote>

<footnote n="15">¹⁵ http://archives.esac.esa.int/hsa/whsa/</footnote>
returned a sample of 33,322 objects. We then queried the Herschel Science Archive with that list to find all observations where the polygon defining the field of view of the observation intersects with a circle of radius of 1' centered at the galaxy's central coordinates. This search was done for all PACS and SPIRE photometric observing modes, aside from parallel modes. The query results in 1245 matches for PACS observations and 1729 for SPIRE (note that some of these are repeated observations of the same galaxies). We also included three very nearby galaxies observed in PACS-SPIRE parallel mode: M31, M33, and IC 342. When a galaxy is observed multiple times, we use all observations (see Section 3.1.2). Some of our targets lie at the edge of Herschel scans and are eliminated from the sample after visual inspection of each image.

Some nearby galaxies were also observed as part of large-area surveys for studying galaxy evolution at higher redshift, as in the large map of the Virgo Cluster (the Herschel Virgo Cluster Survey, HeViCs; J. I. Davies et al. 2010). Due to the computational resources required to reduce large area maps, we do not include those in our reprocessing of archival data. In the case of HeViCs, we extract galaxies from data products already delivered to the Herschel Archive, as described below. Combining the individual pointed observations in PACS or SPIRE; HeViCs; and M31, M33, and IC 342; our initial sample includes 1580 unique galaxies with some Herschel observations.

Although we will only use galaxies with both PACS and SPIRE (as well as the Wide-field Infrared Survey Explorer, WISE) observations in the dust SED fitting described below, we have reduced all galaxies with any PACS or SPIRE maps that are not in a deep field observation and provide those in our data delivery (described in Appendix A).

In the following, we use distances and orientation parameters from the HyperLeda database collated in the galbase code,<fnref n="16" /> updated in some cases to reflect detailed studies of small samples (e.g., using rotation curve fitting). The list of adopted positions, distances, and orientation parameters for all targets in our sample can be found in the delivery table, and a list of the contents of the catalog is provided in Table 1 (adopted parameters are mostly the same as in A. K. Leroy et al. 2019). For the SED-fitting analysis, we narrow our sample to the galaxies with WISE mid-IR observations from the z=0 Multiwavelength Galaxy Synthesis project (z0MGS; A. K. Leroy et al. 2019), as the mid-IR bands are necessary to perform robust fits using a physical dust model. This leads to a data set of 877 galaxies.

In Figure 1 we show our galaxy sample overlaid on the full sample from z0MGS in M\*–SFR space. The z0MGS sample is approximately complete for galaxies above 10⁹ M\*, within 50 Mpc. The top panel of Figure 1 shows the fraction of galaxies in our Herschel sample compared to the z0MGS sample, in bins of stellar mass. That fraction is relatively constant above 10⁹ M\*, where the reference sample is complete. We see that ~10% of the z0MGS sample above 10⁹ M\* is included in our analysis. The SFR–M\* space is filled relatively homogeneously within the completeness bounds z0MGS sample. This suggests that the selection of galaxies with Herschel observations is representative of the population of nearby galaxies, allowing us to draw conclusions about the typical behavior of radiation fields and dust properties for the z ~ 0 galaxy population.

<caption>**Table 1**</caption>

Columns Included in the Delivery Table, Gathering Galaxy Properties Used in This Study, and Integrated Fit Results

<table><thead><tr><th>Parameter</th><th>Unit</th><th>Description</th></tr></thead><tbody><tr><td>GALNAME</td><td>...</td><td>Galaxy “common” identifier</td></tr><tr><td>PGCNAME</td><td>...</td><td>Galaxy PGC identifier, from galbase</td></tr><tr><td>RA_DEG</td><td>deg</td><td>R.A., from galbase</td></tr><tr><td>DEC_DEG</td><td>deg</td><td>Decl., from galbase</td></tr><tr><td>POSANG_DEG</td><td>deg</td><td>Position angle, from galbase; if none, set to 0 during processing</td></tr><tr><td>INCL_DEG</td><td>deg</td><td>Inclination, from galbase; if none, set to 0 during processing</td></tr><tr><td>DIST_MPC</td><td>Mpc</td><td>Distance, from galbase</td></tr><tr><td>R25_DEG</td><td>deg</td><td>Optical radius, from galbase</td></tr><tr><td>REFF_DEG</td><td>deg</td><td>Effective radius, from J. Sun et al. (2024, in preparation)</td></tr><tr><td>FIT_TYPE</td><td>...</td><td>Type of fit: I indicates integrated, and R resolved, which implies integrated fit available as well.</td></tr><tr><td>(e) WISE1</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band WISE 1</td></tr><tr><td>(e) WISE2</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band WISE 2</td></tr><tr><td>(e) WISE3</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band WISE 3</td></tr><tr><td>(e) WISE4</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band WISE 4</td></tr><tr><td>(e) PACS70</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band PACS 70</td></tr><tr><td>(e) PACS100</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band PACS 100</td></tr><tr><td>(e) PACS160</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band PACS 160</td></tr><tr><td>(e) SPIRE250</td><td>Jy</td><td>Integrated flux within r<sub>25</sub> in band SPIRE 250</td></tr><tr><td>(e) UMIN</td><td>...</td><td>Minimum radiation field from the integrated fit</td></tr><tr><td>(e) GAMMA</td><td>...</td><td>Fraction of dust heated by a power-law combination of radiation field from the integrated fit</td></tr><tr><td>(e) UBAR</td><td>...</td><td>Average radiation field, calculated from U<sub>min</sub> and &gamma;</td></tr><tr><td>(e) QPAH</td><td>%</td><td>Fraction of dust mass in the form of PAHs from the integrated fit</td></tr><tr><td>(e) MDUST</td><td>M<sub>&odot;</sub> pc<sup>-2</sup></td><td>Total dust mass from the integrated fit</td></tr></tbody></table>

Note. The (e) denotes the associated error column.

### 3. Data Reduction and Analysis

After the initial sample definition from the query to the Herschel Science Archive, we downloaded the Level 0 data for all bands and all targets and reduced them with the Herschel Interactive Processing Environment (HIPE) and Scanamor-phos pipelines (see Section 3.1). This initial sample includes 1580 galaxies. Note that the SPIRE bands at 350 and 500 µm are not used to derive dust properties (see Section 4), but we have reduced those bands with the same processing steps described below and they are included in our public data delivery. In the following sections, we describe in more detail the data reduction and processing that was applied to all galaxies. We then describe the further data processing needed for the SED-fitting analysis, including the incorporation of the z0MGS WISE data.

#### 3.1. Far-IR: Herschel

##### 3.1.1. Initial Processing

The PACS and SPIRE queries to the Herschel Science Archive were done independently, and we require in both cases that the target was observed in any of the photometry modes. From the queries, we extract the target name and observation

<footnote n="16">¹⁶ https://github.com/akleroy/galbase</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Stellar mass vs. star formation rate for the z0MGS sample (gray points; A. K. Leroy et al. 2019) highlighting the galaxies with WISE and Herschel maps used in this study as black points. Above $10^9 M_\odot$, the z0MGS sample is approximately complete for galaxies within 50 Mpc. The top panel shows the relative number of galaxies in our sample compared to the z0MGS sample in bins of $M_\odot$. We find that ~10% of the z0MGS galaxies have Herschel coverage and are included in our sample. That fraction does not change substantially with stellar mass.</figcaption></figure>

IDs for each source. We processed the parallel-mode observations for M31, M33, and IC 342 in a similar manner, simply changing the relevant keyword in *Scanamorphos* to reduce them as such. We note that for very extended galaxies on the sky, additional processing would be needed to recover very diffuse, extended emission, following C. J. R. Clark et al. (2023). At present, our reduction does not account for any filtered large-scale emission.

We use *HIPE* (version 15; E. Wieprecht et al. 2009; C. D. Dowell et al. 2010; S. Ott 2010) and the *Scanamorphos* suite of routines (version 25, released 2016 August; H. Roussel 2013) to process scans from L0 to L2 products (we use the PACS calibration tree v77, and the SPIRE calibration tree v14.3). This first step takes the raw data, formats it into images, and removes the instrumental noise, as described in H. Roussel (2013). It returns FITS files containing the formatted image (in Jy pix⁻¹ for PACS, and Jy beam⁻¹ for SPIRE) and associated uncertainty, a weight map, and a drift map. The weight map reflects the number of scans per pixel, and informs our assessment of the instrumental noise. We use these maps to weight pixels with respect to this number of scans.

### 3.1.2. Targets with More than One Observation

A few galaxies in our sample were targeted more than once throughout Herschel operations. When this is the case, we manually combine all scans targeting the galaxy into the same formatting framework before running the pipeline. In *Scanamorphos*, this is done by appending scan filenames within the same call to the formatting function. This is valid when observations are done with similar modes. Observations in “mini-maps” mode were included with the appropriate

keywords. We do not include observations that target supernovae to avoid adding in a time variable source.

### 3.1.3. Extracting Galaxies from HeViCS

HeViCS (J. I. Davies et al. 2010; HeVICS Team 2020) is a large program that mapped the Virgo galaxy cluster using the PACS 100 µm through the SPIRE 500 µm bands by observing four contiguous large fields. The reduced HeViCS data have been delivered to HSA as a high-level science product.<fnref n="17" /> HeViCS includes 12 galaxies that have also been observed with targeted observations in other programs. We compare our reprocessing of these 12 galaxies to the HeViCS high-level products by extracting fluxes within $r_{25}$. We find good agreement between the fluxes, which confirms that a full reprocessing of the large HeViCS mosaics is not necessary. As such, we use the products from HeViCS in line with those produced by our reprocessing. For galaxies covered by HeViCS that also have dedicated observations, we use the dedicated observations. We extract a total of 65 galaxies from HeViCS.

To extract galaxies from HeViCS, we select a $5 r_{25}$ area in each band around the central coordinates of each galaxy. In comparing to WISE 4, we found a small astrometric offset between the extracted Herschel and the WISE images for most of these galaxies. We used WISE as a standard, and correct the extracted HeViCS images. This is done by finding the translation (we assume the misalignment is translation only) that maximizes the correlation between the reference image (WISE) and the newly extracted HeViCS cutout. We then update the headers of the images to reflect the adjusted coordinates. After these two steps, we process the HeViCS maps the same way as our other data, including background removal, convolution to a common resolution, and regridding.

## 3.2. Mid-IR: WISE

We combine our far-IR data with the WISE (E. L. Wright et al. 2010) images compiled and analyzed in the z0MGS project (A. K. Leroy et al. 2019). WISE mapped the near- and mid-IR sky at $\lambda \sim 3.4, 4.6, 12$, and $22 \mu m$ (hereafter WISE 1, WISE 2, WISE 3, and WISE 4). We refer to A. K. Leroy et al. (2019) for the detailed description of the catalog creation and image processing. We use their data products at $15''$ resolution, which we will then convolve to match the resolution of the SPIRE 250 µm observations.

In some specific cases, the WISE cutouts provided by A. K. Leroy et al. (2019) are much smaller than the Herschel image and significantly limit the available pixels to measure the background covariance matrix (Section 4.2). For these galaxies, we extract larger cutouts from the WISE tiles to ensure proper measurement of the background. This does not affect any other aspect of the fitting.

## 3.3. Unit Conversion and Background Removal

The output from the initial reduction for PACS and SPIRE observations are in Jy pixel⁻¹ and Jy beam⁻¹, respectively. We convert PACS data to MJy sr⁻¹ using the pixel size. The SPIRE data are converted to MJy sr⁻¹ and corrected from point source to extended source calibration using the $K_{P\text{toE}}$ factors in the SPIRE Handbook v17.<fnref n="18" />

<footnote n="17">¹⁷ https://irsa.ipac.caltech.edu/data/Herschel/HeVICS/overview.html</footnote>

<footnote n="18">¹⁸ http://herschel.esac.esa.int/Docs/SPIRE/spire_handbook.pdf</footnote>
Far-IR observations are subject to several foreground/
background contributions, including MW cirrus, zodiacal light,
and instrumental offsets. To remove these, we perform a
background subtraction with a similar approach to that in
A. K. Leroy et al. (2019) and D. Utomo et al. (2019) using a 2D
plane. For each target, we create a mask that is used to exclude
galaxy emission from the background plane fit. We mask:

1. All pixels within a radius $r = A \times r_{25}$ from the galaxy center. The default value is $A = 1.5$. We require a minimum of 50 pixels to be masked within that radius. If this criteria is not met, $A$ is increased until this inner mask has this minimum required area. Alternatively, $A$ is decreased if the number of unmasked pixels (i.e., background pixels) is less than 10% of the total number of pixels in the image. No bounds are applied to vary $A$, but its final values range between 1 and 1.6, for the whole sample;

2. Other known galaxies in the field of view, identified from HyperLeda, are masked with their own coordinates, and $A \times r_{25}$ as described above;

3. From the weight map given by Level 1 Scanamorphos products, we mask pixels that are below the median value of that weight map. This ensures that we use only the most well-covered pixels to estimate the background;

4. Finally, we perform an iterative clipping, masking any pixel above the median +3× the standard deviation of the background pixels. This is done several times, recomputing median and standard deviation, until the fractional difference between two calculated values of the median, $m_i^{\text{bkg}}$, $m_{i+1}^{\text{bkg}}$, is less than 1%.

The final set of unmasked pixels is fit with a 2D plane
representing the background. This plane is then subtracted from
the image. In our data product delivery, we provide the
background-subtracted image, the final background mask, and
store the coefficients describing the subtracted plane in the
headers.

3.4. *Convolution and Regridding*

We perform convolutions to put all mid- and far-IR images
at matched resolution prior to our pixel-by-pixel SED
modeling. First, all WISE and PACS data are convolved to
the circularized SPIRE 250 point-spread function (PSF;
FWHM ~ 18"). This provides the highest angular resolution
for the dust SED-fitting results, but retains the characteristics of
the SPIRE 250 PSF, which has some non-Gaussianity. For
convolution, we use the circularized kernels from G. Aniano
et al. (2011) for PACS-to-SPIRE convolutions. We also
convolve the native resolution SPIRE 250 to a circularized
SPIRE 250 PSF, following G. Aniano et al. (2011). For WISE,
we use a kernel to convolve a 15" Gaussian to circularized
SPIRE 250 also from G. Aniano et al. (2012). For each target,
we reproject all convolved images to the pixel grid of the
SPIRE 250 image. This grid has pixel size 4"/5 and so
oversamples the SPIRE beam with ~18" FWHM, as has been
done in other works (e.g., D. Paradis et al. 2023). Another
possible approach seen in the literature is to set the pixel size to
that of the PSF. In that case, pixels can be considered
independent, carrying little correlation between them, and
treated as unique data points (e.g., K. D. Gordon et al. 2014;
S. Viaene et al. 2014; N. Vutisalchavakul et al. 2014; V. Casasola
et al. 2015, 2022; G. Saikia et al. 2020). In our case, the analysis is

carried on a large sample that we bin, and extract medians, and the
oversampling will have a limited effect.

3.5. *Background Covariance Matrix Masks*

As described in Section 4, our SED fit makes use of a covariance matrix describing the background at each wavelength included in the fit. To determine this for each galaxy, we need to select which pixels to attribute to the background only, i.e., which are free of any target galaxy or bright foreground star emission (effectively, faint stars and background galaxies are included in the covariance matrix as confusion noise). These pixels are used to measure the background covariance matrix (described in Section 4.2). We construct this mask using the background-subtracted images after convolution and regridding. To identify background pixels in these final images, we follow the same steps used to estimate the background in Section 3.3, adopting the already determined *A* coefficient from that step. The bright pixels identified in any band during iterative masking are excluded from the covariance estimation, as are the regions affected by bright stars in the mid-IR bands in the z0MGS delivery.

3.6. *Integrated Galaxy Photometry*

We use aperture photometry to derive the integrated flux
within $r_{25}$ for each band and each target. We also estimate an
associated uncertainty by scaling down the value of the pixel-
by-pixel noise measured from a signal-free region outside of
$r_{25}$ by the square root of the number of resolution elements in
$r_{25}$. These measurements are included in our delivered data
products, summarized in Table 1.

At this stage, we trim a subset of galaxies from the sample
based on the signal-to-noise of their integrated photometry. If
the integrated flux of any available PACS or SPIRE bands is
below $1\sigma$, we do not perform an integrated SED fitting on that
target (and therefore no resolved fit either). We cut out about
120 targets based on this criterion.

3.7. Validation against Published Far-IR Measurements

In the following, we compare our PACS and SPIRE maps to
previous measurements. Differences may arise from distinct
pipeline processing, background subtraction (particularly at
low surface brightness), and/or extraction apertures. For each
comparison, we quantify the median offset between the fluxes
or surface brightnesses, the $1\sigma$ scatter, and the median absolute
deviation (MAD), all in dex. Figures 2 and 3 present these
comparisons.

KINGFISH. We compare resolved surface brightnesses in our images to those measured by G. Aniano et al. (2020).<fnref n="19" /> They investigated the dust properties in the KINGFISH galaxies (R. C. Kennicutt et al. 2011), using Spitzer and Herschel observations, working at the same SPIRE 250 resolution that we use for our SED fitting. Nearly all KINGFISH galaxies are included in our sample, so the comparison of the resulting maps provides a one-to-one check on our procedures (NGC 5398 and NGC 5408 are not in the GALEX-WISE z0MGS sample, and we omit them here). In Figure 2, we reproject our maps onto the same grid as the G. Aniano et al. (2020) images and compare the intensities between the two images pixel by pixel.

<footnote n="19">¹⁹ The G. Aniano et al. (2020) maps were obtained from http://arks.princeton.edu/ark:/88435/dsp01hx11xj13h.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Resolved comparisons of surface brightnesses, in MJy sr⁻¹, for overlapping galaxies between the KINGFISH (G. Aniano et al. 2020) and Herschel-zOMGS samples. The top row shows the KINGFISH values as a function of the Herschel-zOMGS surface brightness. The 1:1 relation is shown with the black dashed line, and we add offset lines (showing factors of 2 and 5 difference) to guide the eye. The bottom row shows the difference (Herschel-zOMGS-KINGFISH) in dex. We compute offset, scatter, and median absolute deviation (MAD) for the whole sample, shown in the upper-left corner of each panel, and for the “bright” pixels, chosen as log(I<sub>ν</sub><sup>OMGS</sup> [MJy sr<sup>-1</sup>]) ≥ 0.75, shown in the bottom-right corner of each panel.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Comparisons of the integrated photometry from the zOMGS sample compared to LVL (54 galaxies; top row), DustPedia (553 galaxies; bottom row). In general, we find good agreement between the integrated photometry between our measurements and those of DustPedia and LVL, with the largest scatter for fainter galaxies. Average noise values range between 2 and 4 × 10⁻³ Jy in each band, corresponding to ~0.51, 0.36, 0.47, and 0.31 MJy sr⁻¹ from PACS 70 to SPIRE 250.</figcaption></figure>

The maps show outstanding agreement at high intensity, with little scatter, indicating that the processing yields consistent results with KINGFISH. The larger scatter at low intensities (log₁₀(I<sub>ν</sub> [MJy sr<sup>-1</sup>]) ≲ 0.75) is due to the sensitivity of the KINGFISH maps, which is typically 5–7 MJy sr<sup>-1</sup> for PACS and ~1 MJy sr<sup>-1</sup> for SPIRE 250 µm (R. C. Kennicutt et al. 2011).

LVL. D. A. Dale et al. (2023) compiled SEDs, including Herschel photometry for targets of the Local Volume Legacy survey (D. A. Dale et al. 2009). They conducted a careful by-hand definition of local backgrounds and masking of contaminants like stars. In the bottom row of Figure 3, we compare their integrated

fluxes and ours for 54 galaxies that overlap between samples. The overall agreement is excellent, with the integrated fluxes agreeing within a few percent on average and showing typically 0.01–0.02 dex MAD. This overall good agreement with previous measurements confirms the validity of our data processing.

*DustPedia*. We also compare the integrated fluxes measured within *r*<sub>25</sub> for our sample to integrated photometry for the DustPedia sample (C. J. R. Clark et al. 2018), available from the DustPedia archive. An extensive description of the data processing is given in C. J. R. Clark et al. (2018), with details on the CAAPR pipeline (C. J. R. Clark et al. 2015; P. De Vis et al. 2017a, 2017b). In Figure 3, bottom row, we compare
integrated fluxes between our work and DustPedia for 553 galaxies. The agreement between DustPedia and our sample is relatively good, though not as good as our agreement with the LVL and KINGFISH work. We have visually checked the background for galaxies showing a deviation greater than 50% between the two samples in PACS 100, and found no noticeable artifacts. The scatter is likely coming from the different data reduction approaches and differences in the integrated photometry aperture radius. Although the MAD statistics are not as low as with the LVL sample, they are still within a reasonable range, and we consider that the agreement with these two samples is overall satisfactory.

## 3.8. Data from the z0MGS Sample

### 3.8.1. Integrated Values

We compare our observed dust and radiation field properties to galaxy-integrated values of the SFR, stellar mass ($M_*$), and offset from the star-forming main sequence ($\Delta$MS) adopted from the z0MGS catalog (A. K. Leroy et al. 2019). The data were extracted from the z0MGS delivery (z0MGS Team 2019). The definition of the star-forming main sequence used to calculate the offset is given by Equation (19) in A. K. Leroy et al. (2019), reproduced here for clarity:

$$ \log_{10}(\text{sSFR}_{\text{MS}}) [\text{yr}^{-1}] = -0.32 \log_{10} \left( \frac{M_*}{10^{10} M_{\odot}} \right) - 10.17. \quad (1) $$

For each integrated galaxy measurement, we calculate the main-sequence offset as:

$$ \Delta\text{MS} = \log_{10}(\text{sSFR}) - \log_{10}(\text{sSFR}_{\text{MS}}). \quad (2) $$

A. K. Leroy et al. (2019) found that $\log(\text{sSFR}) \le -11$ should not be considered a robust estimate using the combination of UV and IR emission. In the following analysis, we only use targets above this minimum threshold.

### 3.8.2. Galaxy-averaged Surface Densities

To bridge the integrated and resolved measurements, we also calculate “galaxy-averaged” surface densities by dividing the $M_*$ and SFR by the galaxy’s effective area. This calculation uses measurements of effective radii ($r_e$) from J. Sun et al. (2024, in preparation) for the z0MGS sample, including inclination corrections. We divide the integrated SFR and $M_*$ from the A. K. Leroy et al. (2019) atlas by the effective area ($\pi r_e^2$) for each galaxy where an $r_e$ value is available (176 galaxies do not have $r_e$ measurements). We refer to these values as $\langle \Sigma_{SFR} \rangle$ and $\langle \Sigma_* \rangle$.

### 3.8.3. Surface Densities

We compute resolved maps of the SFR surface density ($\Sigma_{SFR}$), stellar mass surface density ($\Sigma_*$), and resolved sSFR, sSFR$_R$, using the prescriptions in A. K. Leroy et al. (2019). To compute the stellar mass surface density, $\Sigma_*$, we use the estimated galaxy-integrated WISE 1 mass-to-light ratio $\Upsilon_*^{3.4}$ for each target provided in the z0MGS catalog and

$$ \Sigma_* \approx 3.3 \times 10^2 \left( \frac{\Upsilon_*^{3.4}}{0.5} \right) I_{3.4 \, \mu\text{m}} \quad (3) $$

where $\Sigma_*$ is in $M_{\odot}$ pc$^{-2}$, and the surface brightness of the WISE 1 map, $I_{3.4 \, \mu\text{m}}$, is in MJy sr$^{-1}$. While the true appropriate

mass-to-light conversion may vary from location to location within a galaxy, the choice to use a single galaxy-integrated value ensures that the integrated galaxy results and the resolved results yield consistent stellar masses.

To compute $\Sigma_{SFR}$, we use the WISE 4 photometry combined with GALEX FUV or NUV (with the top priority being FUV, next NUV if FUV is not available, and WISE 4 alone otherwise), per the hybrid prescription in A. K. Leroy et al. (2019):

$$ \Sigma_{SFR} \approx f_{\text{FUV,NUV}} \left( \frac{C_{\text{FUV,NUV}}}{10^{-43.35}} \right) I_{\text{FUV,NUV}} + 3.24 \times 10^3 \left( \frac{C_{\text{WISE 4}}}{10^{-42.7}} \right) I_{22 \, \mu\text{m}}, \quad (4) $$

with $I_{FUV,NUV}$, $I_{22 \, \mu\text{m}}$ in MJy sr$^{-1}$, $f_{FUV} = 1.04 \times 10^{-1}$, and $f_{NUV} = 1.05 \times 10^{-1}$. We refer to Table 7 in A. K. Leroy et al. (2019) to find the C coefficients to use to convert the GALEX FUV or NUV and WISE 4 maps to $\Sigma_{SFR}$ in $M_{\odot}$ kpc$^{-2}$ yr$^{-1}$.

We calculate the resolved sSFR as $sSFR_R = \Sigma_{SFR}/\Sigma_*$. To compute the resolved maps of the offset from the main sequence ($\Delta$MS$_R$), we first define the resolved star-forming main sequence using the maps of $\Sigma_*$ and sSFR$_R$. We fit a line to binned averages of all pixels above the completeness thresholds (see Section 4.4). We find $\log_{10}(\text{sSFR}_{MS,R}) = -0.625 \times \log_{10}(\Sigma_*) - 5.216$. We use that resolved main sequence to then calculate the offset from that prediction of all pixels based on their resolved sSFR, as

$$ \Delta\text{MS}_R = \log_{10}(\text{sSFR}_R) - \log_{10}(\text{sSFR}_{MS,R}). \quad (5) $$

We note that there are a variety of methodological choices possible for defining a resolved main sequence, so Equation (5) should be viewed primarily as a characterization of the offset from the average sSFR$_R$ at a given $\Sigma_*$ specifically for our sample.

## 4. Dust Emission Modeling

We fit the WISE through Herschel 250 µm SEDs of dust emission using the B. T. Draine & A. Li (2007) physical dust models for both resolved and integrated measurements of each galaxy in our sample. We work at SPIRE 250 µm resolution because it provides sufficient coverage of the long wavelength dust SED, while also preserving angular resolution (G. Aniano et al. 2012; J. Chastenet et al. 2021). In Appendix C, we show the dust parameter distributions from fits using IR data up to SPIRE 500 (x-axes) against those from fits up to SPIRE 250 (y-axes), from the work of J. Chastenet et al. (2021). This shows that trends and parameter values are well kept without using the last two SPIRE bands.

### 4.1. The B. T. Draine & A. Li (2007) Model

We use the B. T. Draine & A. Li (2007) model with the dust mass renormalization factor of 3.1 determined by J. Chastenet et al. (2021), which results from fitting the Milky Way high-latitude cirrus spectrum from M. Compiègne et al. (2011) with depletion constraints from E. B. Jenkins (2009) and K. D. Gordon et al. (2014). We use the B. T. Draine & A. Li (2007) MW (“Milky Way”) model, with $R_V = 3.1$. Our procedure follows that of J. Chastenet et al. (2021) exactly, and all details in this section can also be found there in more detail.
The radiation field heating the dust is described by a model with a delta function at a minimum radiation field, $U_{min}$, and a power-law distribution of radiation field intensities extending to $U_{max}$ (D. A. Dale et al. 2001). The radiation field intensity is scaled in units of the J. S. Mathis et al. (1983) solar neighborhood radiation field intensity, $U_{MMP}$. The delta function plus power-law model has been used in a variety of dust SED modeling approaches, particularly with the B. T. Draine & A. Li (2007) model (e.g., B. T. Draine et al. 2007; G. Aniano et al. 2012, 2020). In this model, a fraction $(1-\gamma)$ of the dust grains is illuminated by $U_{min}$, while the fraction $\gamma$ is heated by the power-law distribution with $U_{min} < U \le U_{max}$, defined as

$$ \frac{1}{M_{\mathrm{d,tot}}} \left( \frac{dM_{\mathrm{d}}}{dU} \right) = (1-\gamma)\delta(U - U_{\mathrm{min}}) + \gamma \frac{\alpha - 1}{U_{\mathrm{min}}^{1-\alpha} - U_{\mathrm{max}}^{1-\alpha}} U^{-\alpha}. \quad (6) $$

We fix the power-law coefficient to $\alpha=2$, and we use $U_{max} = 10^7$ (following G. Aniano et al. 2012, 2020). Using the fit results for $U_{min}$ and $\gamma$, we also compute the average radiation field, $\bar{U}$:

$$ \bar{U} = (1 - \gamma) U_{\mathrm{min}} + \gamma \times U_{\mathrm{min}} \frac{\ln(U_{\mathrm{max}}/U_{\mathrm{min}})}{1 - U_{\mathrm{min}}/U_{\mathrm{max}}}. \quad (7) $$

We fit for the total dust surface density, $\Sigma_d$, which combines both carbonaceous and silicate grain populations, and the mass fraction of carbonaceous grains with $<10^3$ carbon atoms, $q_{PAH}$. We model the mid-IR starlight contribution to the SED (most notable at WISE 1 and 2) with a scaled 5000 K blackbody, $\Omega_*$.

## 4.2. Fitting Procedure

We use the grid-based Bayesian fitting code *DustBFF* (K. D. Gordon et al. 2014) to fit the dust emission model to the data. We use flat priors for all five parameters described in the previous section. Table 2 lists the fitted parameters, with the sampling range and step size of the grid.

*DustBFF* uses a multivariate Gaussian with covariance matrix $C$ to model the uncertainties due to statistical noise, background flux, and instrument calibration and their correlations. We use the same values as J. Chastenet et al. (2021) for the instrumental uncertainties. The background covariance matrix measures the correlation of the background signal for each band (Equation (23) in K. D. Gordon et al. 2014). It changes for each set of mid- through far-IR images, and therefore we build a unique matrix for each galaxy (Section 3.5). For fits to the integrated SEDs, the integrated fluxes would not be appropriately modeled with a background covariance matrix calculated with the same approach. In addition, the instrumental calibration uncertainties dominate over noise for the majority of galaxies. We therefore omit the background covariance matrix in this case.

The fitting is done at the SPIRE 250 resolution (~18"). We fit all pixels that are above $1\sigma$ of the background in each band (not including WISE 1 and WISE 2), for a greater spatial coverage, and possible future use. However, we only use pixels above $3\sigma$ of the background in each band (not including WISE 1 and WISE 2) for the following figures and analysis, as they offer more reliable results.

<table><caption>**Table 2**
Fitting Parameters</caption><thead><tr><th>Parameter</th><th>Range</th><th>Step</th><th>Unit</th></tr></thead><tbody><tr><td>log<sub>10</sub>(Σ<sub>d</sub>)</td><td>[-3.0, 0.7]</td><td>0.1</td><td>M<sub>☉</sub> pc<sup>-2</sup></td></tr><tr><td>q<sub>PAH</sub></td><td>[0, 6.1]</td><td>0.25</td><td>%</td></tr><tr><td>U<sub>min</sub></td><td>[0.1, 50]</td><td>Irr<sup>a</sup></td><td>...</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>[-4.0, 0.0]</td><td>0.15</td><td>...</td></tr><tr><td>log<sub>10</sub>(Ω<sub>*</sub>)</td><td>[-1.5, 3.0]</td><td>0.1</td><td>...</td></tr></tbody></table>

**Note.**

a $U_{min} \in \{0.1, 0.12, 0.15, 0.17, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.2, 1.5, 1.7, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 12.0, 15.0, 17.0, 20.0, 25.0, 30.0, 35.0, 40.0, 50.0\}$.

## 4.3. Integrated and Resolved Fits

We perform dust emission fitting on both integrated and resolved measurements. Our final samples include a total of 819 targets with resolved fits. We also fit the integrated SED for an additional 58 targets that were not suited for pixel-by-pixel fitting due to their low surface brightness but were still significantly detected. This brings the number of galaxies with integrated fits to 877.

## 4.4. Pixels Used in the Analysis

We performed fitting on all pixels that passed a $1\sigma$ cut in all of the relevant infrared bands (~734,800 total pixels corresponding to ~41,000 independent resolution elements). For further analysis, we make a series of stricter cuts to select high-confidence dust emission measurements. We also make several quality cuts to remove data with artifacts and poor fits from the sample. These cuts include:

1. *S/N > 3 in infrared bands.* Although the fit is performed on all pixels detected at $1\sigma$ above the background in all bands from WISE 3 to SPIRE 250, we only use pixels passing a $3\sigma$ cut in all bands in the following analysis. This ensures reasonable confidence in the SED-fitting results for each pixel. Out of the original sample, the number of pixels that pass this cut is ~474,900 corresponding to ~26,000 independent resolution elements.

2. *Saturation in WISE.* For bright nuclear point sources and for some very high surface brightness starburst galaxies, the WISE observations can be saturated. The z0MGS catalog provides information on saturation in the WISE bands for their whole sample. We use their surface brightness thresholds (at the 7"/5 resolution), 100 MJy sr$^{-1}$ in WISE 1 and WISE 2 and 300 MJy sr$^{-1}$ in WISE 3 and WISE 4, to cut out pixels exceeding these values, which are likely to experience saturation. We applied these “saturation masks” to NGC 598, NGC 1097, NGC 1365, NGC 3621, NGC 5236, and NGC 6946. Upon visual inspection, all other galaxies flagged with saturation in the z0MGS catalog were discarded from the analysis due to having more extensive artifacts related to saturation. For IC 342, we also perform by-hand masking, discarding pixels that show unrealistically high sSFR$_R > 10^{-8.5}$ yr$^{-1}$, likely due to foreground star contamination. These cuts remove a total of ~30,000 pixels (2000 independent data points).

3. *SED fit residuals.* We only keep pixels that show fit residuals (described as $(\text{data} - \text{model})/\text{model}$) below 50% in all PACS bands. This choice is discussed
in further detail below in Section 4.5. This cuts out 240 pixels.

Unless mentioned otherwise, all subsequent figures only show pixels that pass these quality cuts.

## 4.5. Radiation Field Parameters and Limits of the Fitting

The B. T. Draine & A. Li (2007) models span a range of possible radiation field properties, with $U_{min}$ ranging between 0.1 and 50 times the $U_{MMP}$ radiation field and a potential range of $log(\gamma)$ from -4 to 0. While this range encompasses the vast majority of points we expect to encounter in nearby star-forming galaxies, there are some environments where the radiation field may not be well represented by our models.

In intense starbursts, where the radiation field on ~kpc scales can be >50 $U_{MMP}$, the models cannot reproduce the pervasive, high, minimum radiation field intensity. In such cases, we find that the fit results suggest a best fit where more of the dust luminosity is shifted into the power-law distribution. In this case, we typically see an underprediction of the 70–100 µm emission, due to the inability of the models to extend to $\bar{U}$ values that generate SEDs peaking at 70–100 µm. Such issues manifest as negative residuals from the fit in the PACS 70 and 100 bands, and inspection shows they are correlated with high fit $\gamma$ values.

At the other end of the radiation field intensity range, it is important to note that we only include observations out to 250 µm in our fit. This means that for low radiation field intensities, where the peak of the thermal dust emission moves close to or beyond 250 µm, we have a poor handle on the dust temperature and therefore dust mass, leading to bad fits. Selecting galaxies with $\langle\bar{U}\rangle \sim \langle U_{min} \rangle \le 0.3$ (corresponding to $T \sim 15$ K), at the low end of our range of available $U_{min}$, all show large residuals in PACS 100, PACS 160, and SPIRE 250. This suggests that these galaxies have dust cold enough that the fit up to 250 µm cannot reproduce the SEDs with the available $U_{min}$ in our models (lower limit of $U_{min}$ is 0.1). Additionally, we note that some of these galaxies do not have PACS 70 data.

We find that both of these failure cases can be flagged via a cut on residuals in the PACS bands. We find that a cut on the residuals in the PACS bands at $\text{abs}((\text{data} - \text{model})/\text{model}) > 0.5$ is effective at isolating fits that show either very high or very low radiation field intensities with poor fits. This cut has a minimal impact on the overall sample, flagging 240 pixels overall, and we have verified that the main effect is to eliminate outliers and none of the scaling relationships change significantly if we tighten this cut to 25% or relax the cut entirely.

## 4.6. Completeness

One of the goals of the analysis we present below is determining scaling relationships between $\Sigma_*$, $\Sigma_{SFR}$, $\Sigma_d$, $q_{PAH}$, and the radiation field properties in our galaxies. In order to accurately measure these relationships, we first need to ensure that our sample is not biased due to the detection limits for Herschel and WISE, for the resolved scaling relations. Because dust emission generally gets fainter as a function of galactocentric radius, our $3\sigma$ cut translates into a cut in $\Sigma_*$ and $\Sigma_{SFR}$ as well. This means that at lower values of $\Sigma_*$ and $\Sigma_{SFR}$, we will be biased toward selecting positive outliers in dust emission, either regions with more dust for a given combination of $\Sigma_*$ or $\Sigma_{SFR}$, or regions where noise pushes points above $S/N = 3$. This results in an overall bias in what we judge to be the

average properties of the dust and radiation field. To quantify and minimize that bias, we calculate completeness thresholds in $\Sigma_{SFR}$, $\Sigma_*$, $sSFR_R$, and $\Delta MS_R$. To do this, we use the distribution of all of the pixels in the resolved maps prior to any cuts, and bin them by $\Sigma_*$, $\Sigma_{SFR}$, $sSFR_R$, and $\Delta MS_R$. In each bin, we calculate the fractions of pixels that pass the $3\sigma$ cut in bands WISE 3–SPIRE 250 (note that the other cuts for saturation and residuals affect only a small fraction of the pixels and so do not matter in judging completeness). We then find the values of $\Sigma_*$, $\Sigma_{SFR}$, $sSFR_R$, and $\Delta MS_R$ where 50% of pixels in a bin are below this $3\sigma$ cut in any of those bands; we define this bin as our completeness threshold.

From the results of this analysis, we consider coverage of the sample complete if

1. $\log_{10}(\Sigma_{SFR} [M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}]) > -2.57$,

2. $\log_{10}(\Sigma_* [M_{\odot} \text{pc}^{-2}]) > 1.82$,

3. $\log_{10}(sSFR_R [\text{yr}^{-1}]) > -10.39$,

4. $\Delta MS_R [\text{dex}] > -0.29$.

These thresholds are marked with vertical white-dashed lines in resolved-fits figures (showing quantities as a function of $\Sigma_*$, $\Sigma_{SFR}$, $sSFR_R$, or $\Delta MS_R$). The analysis and subsequent scaling relations are inferred using only the pixels above these completeness thresholds. This corresponds to ~284,000, 180,000, 285,000, and 295,000 pixels above the thresholds in $\Sigma_{SFR}$, $\Sigma_*$, $sSFR_R$, and $\Delta MS_R$, respectively.

## 4.7. Resolved Fits: Comparison with KINGFISH

In Figure 4, we compare our resolved results to those from G. Aniano et al. (2020) obtained fitting the B. T. Draine et al. (2007) models to the KINGFISH sample, also working at SPIRE 250 µm resolution.

The dust surface densities ($\Sigma_d$) agree well between the two data sets, with a slight nonlinearity in the slope. This likely reflects the adjusted dust opacities that we adopt. In this work, we adopt dust mass renormalization factors for the B. T. Draine & A. Li (2007) models that were derived in J. Chastenet et al. (2021). Compared to the original B. T. Draine & A. Li (2007) model, the renormalization resulted in a factor of 3.1 lower dust mass with no dependence on $U_{min}$ or other radiation field parameters. In G. Aniano et al. (2020), a different set of corrections were applied to the B. T. Draine & A. Li (2007) model based on $\Sigma_d$ offsets measured by comparing all-sky WISE and Planck SED modeling to extinction measured toward background quasars (Planck Collaboration et al. 2016). The correction factor derived in that work included a $U_{min}$ dependence, and the median value over the KINGFISH sample results in 1.6 times lower dust mass. Thus, we might expect to find a factor of approximately 3.1/1.6 ~ 2 lower dust masses in our results based on these different opacity assumptions, with closer agreement at the $U_{min}$ values less than 1, where the G. Aniano et al. (2020) correction factor increases. This generally agrees with what we observe in the comparison of the $\Sigma_d$ results in the left panel of Figure 4, where our results lie below the one-to-one line at the highest $\Sigma_d$ (which tend to occur in inner regions of galaxies where $U_{min}$ is higher) and get closer to the one-to-one line as $\Sigma_d$ decreases. Thus, we argue that our results are in good agreement with G. Aniano et al. (2020) given the different opacity assumptions.

In the two right panels of Figure 4, we show the comparison between $U_{min}$ and $\bar{U}$ between our results and KINGFISH. Here we also see good agreement. Both $\bar{U}$ and $U_{min}$ are offset to
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Comparison of dust and radiation field parameters between Herschel-z0MGS and G. Aniano et al. (2020) for the KINGFISH galaxies (without NGC 5398 and NGC 5408). There is a good agreement between the two samples, despite different fitting approaches. The parameters and units are shown in each panel, with the offset, scatter, and MAD measured with respect to the 1:1 line (black dashed line), except for $q_{\text{PAH}}$ (second panel), where we find a best-fit slope of 0.49, shown in gray. Note that the $\Sigma_d$ distribution crosses the 1:1 line, which leads to a different offset if we measure it below or above a threshold of $\log_{10}(\Sigma_d^{\text{70MGS}}) = -0.6$.</figcaption></figure>

0.1–0.2 dex lower values in our fits relative to KINGFISH. The fit in G. Aniano et al. (2020) allows the slope of the radiation field power-law distribution, $\alpha$, to vary, while we fix it to $\alpha = 2$. While they find that on average $\alpha$ is consistent with ~2, they do find small systematic variations between the inner and outer disks of some galaxies, with $\alpha$ being lower in the centers. The fact that our $\bar{U}$ and $U_{\min}$ are in good agreement suggests that allowing $\alpha$ to vary does not have a significant impact on the fitting and allows us to more efficiently explore the grid space of other fit parameters.

The largest differences in comparison to KINGFISH occur in $q_{\text{PAH}}$. We find that our results yield higher values of $q_{\text{PAH}}$ by ~1% at the highest $q_{\text{PAH}}$. At lower values of $q_{\text{PAH}}$, our results are in better agreement with KINGFISH. These differences likely arise from the fact that we use mid-IR constraints from WISE rather than Spitzer. In general, the WISE 1 and 2 bands (i.e., 3.6 and 4.5 µm) are very comparable to the Spitzer/IRAC 1 and 2 bands, and the Spitzer/MIPS 24 µm is a good match for WISE 22 µm (WISE 4). The primary difference in replacing Spitzer with WISE is therefore the shift of the main constraint on PAHs from the 7.7 µm complex (sampled by the Spitzer/IRAC 8 µm filter) to the 11–12 µm complex (sampled by the WISE 12 µm filter). The WISE 12 µm filter also samples a wider wavelength range than Spitzer/IRAC 8 µm and therefore includes more continuum. The B. T. Draine & A. Li (2007) model includes a fixed PAH spectrum with no variation in the intrinsic band ratios. It may be the case that the difference in $q_{\text{PAH}}$ when using IRAC versus WISE could arise from our targets having brighter 11.3 µm features relative to 7.7 µm compared to what is in the model. In addition, many SINGS/KINGFISH galaxies show evidence for IRAC scan-direction artifacts in the $q_{\text{PAH}}$ maps from G. Aniano et al. (2020), which may reflect some persistence or saturation effects over the bright centers of galaxies. The WISE maps are better behaved in this regard. In addition, the all-sky coverage of WISE also allows for an improved estimation of the background. The Spitzer/IRAC images in particular have relatively narrow fields of view, which could, in some cases, lead to oversubtraction of the diffuse emission in the mid-IR, and underestimation of $q_{\text{PAH}}$. In general, the good correlation of the $q_{\text{PAH}}$ results from the two studies, despite their offset, give us confidence that our fits are tracing $q_{\text{PAH}}$. Future work that incorporates both IRAC and WISE into the fitting or makes use of the more finely sampled mid-IR filter set from JWST may be able to reveal how PAH properties influence $q_{\text{PAH}}$.

To summarize, aside from a few minor differences and an overall scaling of $q_{\text{PAH}}$, we find good agreement between our fitting results and those of G. Aniano et al. (2020). We proceed in the following sections to discuss the results from our SED fitting.

# 5. Results

## 5.1. Radiation Field Scaling Relations: Integrated Fits

We show the distribution of integrated dust and radiation field parameters $U_{\min}$, $\gamma$, and $\bar{U}$ as functions of SFR, $M_*$, sSFR, and $\Delta$MS in Figure 5, and similar scatter plots for $q_{\text{PAH}}$ and $M_d$ in Figure 6. We show all fitted galaxies in light gray, and the darker points are galaxies with $\log_{10}(\text{sSFR}) > -11$. In each panel, we show the binned medians considering only the dark points, with error bars showing one standard deviation. We estimate scaling relations by fitting the medians, where there are at least 10 points in a bin. We quote the results of power-law fits in Table 3.<fnref n="20" /> We note that in general, $U_{\min}$ and $\bar{U}$ behave very similarly, indicating that the average interstellar radiation field is dominated by the $U_{\min}$ component, not the power-law distribution.

These plots reveal several trends:

1. All radiation field parameters ($U_{\min}$, $\bar{U}$, $\gamma$) are positively correlated with quantities involving star formation, i.e., SFR, sSFR, and $\Delta$MS.

2. The correlations with sSFR show the steepest slopes for all radiation field parameters. $U_{\min}$ and $\bar{U}$ both change by around 1 order of magnitude or more as sSFR increases from $10^{-11}$ to $10^{-9}$ yr$^{-1}$.

3. Because the “main sequence” is not flat in $M_*$–sSFR space, it is interesting to check for trends perpendicular to the main sequence at a given $M_*$, quantified by $\Delta$MS. In the right panels of Figure 5, we show that the trends of radiation field parameters $U_{\min}$ and $\bar{U}$ with $\Delta$MS show slightly shallower slopes than the correlation with sSFR.

<footnote n="20">²⁰ Throughout the paper we report power-law fit parameters for each scaling relation. We note that there is no fundamental reason to expect power-law behavior between these variables, and in several cases there is evidence for variations that are not captured by our fits. As an example, the $\bar{U}$–sSFR correlation in Figure 5 shows a relatively flat distribution at $\log(\text{sSFR}) \sim -10$ and a steeper slope above. We argue that to first order the power law provides a straightforward and reasonably accurate representation of the data.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Distribution of $U_{\min}$ (top row), $\gamma$ (middle row), and $\bar{U}$ (bottom row) from modeling integrated galaxies, as a function of sSFR (first column), SFR (second column), $M_*$ (third column), and $\Delta$MS (fourth column). Shown are all of the integrated fits with acceptable residuals (which includes the light-gray points), and the darker points show those with $\log_{10}(\text{sSFR}) > -11$. The white circles are the running medians (on dark points only), with error bars as one standard deviation. The filled circles are those used to fit the scaling relations. The boxed text shows $a \pm \epsilon_a^{\text{fit}} \pm \epsilon_a^{\text{bin}}$, with $Y = a \times X + b$, where $\epsilon_a^{\text{fit}}$ is the error on the slope for the shown fit, and $\epsilon_a^{\text{bin}}$ is the standard deviation of slopes fit with a varying number of bins. These values are tabulated in Table 3.</figcaption></figure>

4. Correlations of $U_{\min}$, $\bar{U}$, and $\gamma$ with SFR all show shallower slopes compared to trends with sSFR and $\Delta$MS.

5. $U_{\min}$, $\bar{U}$, and $\gamma$ have little correlation with $M_*$. This is particularly notable since $M_*$ spans >3 orders of magnitude over the completeness range of our sample.

6. Judged by the root mean square error (RMSE) of the correlation, $\bar{U}$, $U_{\min}$, and $\gamma$ can be best predicted using $\Delta$MS, though SFR also predicts $\gamma$ with nearly the same RMSE. Note that calculating $\Delta$MS requires knowing $M_*$ and SFR.

In general, these findings agree with expectations for how the “delta-function plus power-law model” should behave, and give insights into what is setting the $U_{\min}$ value, which is the dominant part of the ISRF. In theory, both recent star formation and existing older stellar populations (B stars and later) should contribute to setting $U_{\min}$. The power-law distribution is expected to arise from dust heated by higher-intensity radiation fields in the vicinity of recently formed stars. The correlation of $\gamma$ with SFR, sSFR, and $\Delta$MS agrees with this picture. Interestingly, $U_{\min}$ (and $\bar{U}$)

correlations show little dependence on $M_*$, which may indicate that current star formation dominates the radiation field on integrated scales. However, the correlation of $\bar{U}$ and $U_{\min}$ with SFR is also fairly weak, which could be interpreted as extensive quantities like $M_*$ and SFR not capturing the physics that set the radiation field (e.g., the radiation field distribution should differ between a compact, high-intensity distribution of SF compared to a larger, low-intensity distribution). These integrated relationships may also simply reflect a different spatial distribution of the stellar mass, star formation, and ISM throughout the galaxy, for example, if stellar mass is highly centrally concentrated while dust and star formation are more distributed through the disk. We address this topic in Sections 5.2 and 5.5 by looking at the resolved and galaxy-averaged versions of these plots.

## 5.2. Radiation Field Scaling Relations: Resolved Fits

In Figure 7, we show the 2D density histograms of dust parameters $U_{\min}$, $\gamma$, and $\bar{U}$ as functions of sSFR$_R$, $\Sigma_{\text{sSFR}}$, $\Sigma_*$, and $\Delta$MS$_R$, and similar histograms for $q_{\text{PAH}}$ and $\Sigma_d$ in Figure 8. In
<figure><img src="image_10.png" /><figcaption>Figure 6. Same as Figure 5, for fitted parameters $q_{PAH}$ (top row) and $M_d$ (bottom row).</figcaption></figure>

each panel, we show the binned medians, with error bars showing one standard deviation. The filled circles mark the medians used to fit the power-law scaling relationships, selected to have at least 250 measurements in a bin. We report the fits in Table 4 and also show a comparison between the slopes derived for the scaling relationships in the integrated and resolved fits in Figure 9. We note the excellent agreement between the completeness thresholds defined in Section 4.4, computed independently, and the behavior of the binned medians, which begin to noticeably deviate from the trends below the completeness limits. In Appendix B we present correlations between the resolved fit parameters ($\Sigma_d$, $q_{PAH}$, $U_{min}$, $\bar{U}$, and $\gamma$) to give further insight into their behavior.

The resolved scaling relations between radiation field parameters ($\bar{U}$, $U_{min}$, and $\gamma$) and $\Sigma_*$, $\Sigma_{SFR}$, $\text{sSFR}_R$, and $\Delta\text{MS}_R$ reveal the following trends:

1. All radiation field parameters are positively correlated with quantities related to star formation: $\Sigma_{SFR}$, $\text{sSFR}_R$, and $\Delta\text{MS}_R$, similar to what was found for integrated measurements.

2. The steepest slopes occur in the relationships between the radiation field quantities and $\text{sSFR}_R$ or $\Delta\text{MS}_R$, suggesting sSFR and related quantities are important for setting the radiation field distribution.

3. Despite the lack of correlation with $M_*$ in the integrated measurements, $\bar{U}$ and $U_{min}$ do show a significant correlation with $\Sigma_*$ in the resolved measurements. This could result from the older stellar distribution contributing to setting $U_{min}$, with $\Sigma_*$ tracing the density of old stars, or it might reflect that $\Sigma_*$ and $\Sigma_{SFR}$ are often both radially decreasing and therefore correlated.

4. The fraction of dust luminosity from the power-law distribution, $\gamma$, does not strongly depend on $\Sigma_*$, a weakly negative correlation. This reinforces that the older stellar

mass distribution is not critical in setting the power-law part of the radiation field distribution.

5. In general, radiation field parameters are slightly better predicted by $\Delta\text{MS}_R$ compared to $\text{sSFR}_R$. There is a small decrease in the RMSE for all relationships of radiation field quantities with $\Delta\text{MS}_R$ compared to $\text{sSFR}_R$, suggesting the perpendicular offset from the resolved main sequence is better at predicting the radiation field properties than $\text{sSFR}_R$ alone. Given that our measurement of the resolved MS has a nonlinear slope of -0.625 (see Equation (5)), $\text{sSFR}_R$ and $\Delta\text{MS}_R$ are not simply scaled versions of each other, so the slight but consistent preference for $\Delta\text{MS}_R$ as a radiation field predictor suggests that this is meaningful.

6. Judged by the RMSE of the correlation, $\bar{U}$ and $U_{min}$ can be best predicted using $\Sigma_{SFR}$ or $\Delta\text{MS}_R$. $\gamma$ is best predicted with $\Delta\text{MS}_R$. We note that some part of the RMSE can be due to behavior not captured by our power-law fits.

The difference in slope between the resolved and integrated fits potentially provides some insight into the drivers of the radiation field distribution in galaxies. Figure 9 summarizes the slopes of the various scaling relationships as a function of integrated or resolved stellar mass, SFR, sSFR, and $\Delta\text{MS}$. Some trends we see in this comparison include:

1. The slopes of scaling relations for all radiation field quantities ($U_{min}$, $\bar{U}$, and $\gamma$) versus sSFR or $\text{sSFR}_R$ remain consistent within their uncertainties between the integrated to resolved analysis. This may be expected since sSFR is a ratio of surface densities, so may to first order remove radial variation of $\Sigma_*$ and $\Sigma_{SFR}$ within galaxies that leads to differences in resolved and integrated fits.

2. The slopes of $\gamma$ versus $\Sigma_*$ and $M_*$ are slightly different, though both similar to a flat or slightly decreasing trend,
<table><caption>**Table 3**
Coefficients to the Integrated Galaxy Fits, $Y = aX + b$, Presented in Figures 5 and 6</caption><thead><tr><th rowspan="2">Y = Dust Parameter</th><th colspan="5">Binned Medians</th><th colspan="2">All Points</th></tr><tr><th>a ± ε<sub>a</sub><sup>fit</sup> ± ε<sub>a</sub><sup>bin</sup></th><th>b ± ε<sub>b</sub><sup>fit</sup> ± ε<sub>b</sub><sup>bin</sup></th><th>Pearson's ρ</th><th>RMSE</th><th>p-value</th><th>a'</th><th>p-value'</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(SFR)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.21 ± 0.04 ± 0.03</td><td>0.4 ± 0.04 ± 0.03</td><td>0.88</td><td>0.437</td><td>0.002</td><td>0.17</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.27 ± 0.01 ± 0.03</td><td>-1.89 ± 0.02 ± 0.01</td><td>0.99</td><td>0.543</td><td>&lt;0.001</td><td>0.37</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.25 ± 0.04 ± 0.02</td><td>0.53 ± 0.05 ± 0.02</td><td>0.88</td><td>0.444</td><td>0.002</td><td>0.25</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-0.72 ± 0.29 ± 0.24</td><td>2.42 ± 0.36 ± 0.22</td><td>-0.28</td><td>1.898</td><td>0.5</td><td>-0.03</td><td>0.72</td></tr><tr><td>log<sub>10</sub>(M<sub>d</sub> [M<sub>☉</sub>])</td><td>0.81 ± 0.06 ± 0.04</td><td>7.25 ± 0.06 ± 0.02</td><td>0.98</td><td>0.551</td><td>&lt;0.001</td><td>0.81</td><td>&lt;0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(M<sub>*</sub>)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>-0.13 ± 0.06 ± 0.03</td><td>1.56 ± 0.57 ± 0.32</td><td>-0.62</td><td>0.46</td><td>&lt;0.001</td><td>-0.02</td><td>0.58</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.12 ± 0.03 ± 0.04</td><td>-3.08 ± 0.32 ± 0.45</td><td>0.73</td><td>0.581</td><td>&lt;0.001</td><td>0.26</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>-0.06 ± 0.05 ± 0.03</td><td>0.97 ± 0.51 ± 0.33</td><td>-0.4</td><td>0.492</td><td>0.2</td><td>0.05</td><td>0.1</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>1.19 ± 0.26 ± 0.07</td><td>-8.46 ± 2.48 ± 0.71</td><td>0.88</td><td>1.644</td><td>&lt;0.001</td><td>0.77</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(M<sub>d</sub> [M<sub>☉</sub>])</td><td>1.09 ± 0.05 ± 0.04</td><td>-3.52 ± 0.53 ± 0.38</td><td>0.99</td><td>0.466</td><td>&lt;0.001</td><td>1.0</td><td>&lt;0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(sSFR)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.54 ± 0.05 ± 0.03</td><td>5.74 ± 0.54 ± 0.28</td><td>0.95</td><td>0.394</td><td>&lt;0.001</td><td>0.06</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.47 ± 0.03 ± 0.03</td><td>2.77 ± 0.35 ± 0.3</td><td>0.97</td><td>0.563</td><td>&lt;0.001</td><td>0.1</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.7 ± 0.06 ± 0.06</td><td>7.51 ± 0.58 ± 0.64</td><td>0.93</td><td>0.394</td><td>&lt;0.001</td><td>0.09</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-2.58 ± 0.23 ± 0.14</td><td>-22.86 ± 2.28 ± 1.4</td><td>-0.92</td><td>1.541</td><td>&lt;0.001</td><td>-0.32</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(M<sub>d</sub> [M<sub>☉</sub>])</td><td>0.08 ± 0.1 ± 0.04</td><td>7.94 ± 0.97 ± 0.42</td><td>0.18</td><td>0.836</td><td>0.5</td><td>-0.11</td><td>0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = ΔMS</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.44 ± 0.05 ± 0.05</td><td>0.3 ± 0.03 ± 0.02</td><td>0.93</td><td>0.372</td><td>&lt;0.001</td><td>0.46</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.51 ± 0.04 ± 0.04</td><td>-1.95 ± 0.02 ± 0.01</td><td>0.96</td><td>0.522</td><td>&lt;0.001</td><td>0.6</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.57 ± 0.06 ± 0.06</td><td>0.41 ± 0.04 ± 0.03</td><td>0.94</td><td>0.352</td><td>&lt;0.001</td><td>0.6</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-1.89 ± 0.26 ± 0.28</td><td>3.2 ± 0.16 ± 0.08</td><td>-0.89</td><td>1.578</td><td>&lt;0.001</td><td>-1.25</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(M<sub>d</sub> [M<sub>☉</sub>])</td><td>0.67 ± 0.08 ± 0.05</td><td>7.26 ± 0.05 ± 0.04</td><td>0.92</td><td>0.752</td><td>&lt;0.001</td><td>0.59</td><td>&lt;0.001</td></tr></tbody></table>

**Note.** The quoted (a, b) values are for the displayed number of bins, which are running medians of the whole sample; the first uncertainty, $ε_a^{fit}$, is the statistical error on these linear fits; the second uncertainty, $ε_a^{bin}$, is the standard deviation on coefficients (a, b) when fitting binned data with a varying number of bins. Only the filled medians are considered for the fit, and chosen to have at least 10 points within the bin. The Pearson's $\rho$ coefficient is that of the fits shown. The root mean square error (RMSE) is calculated on all dark-gray points. The last two columns show $a'$ and its $p$-value, the slope derived by fitting *all points*, without binning, and give an additional sense of the uncertainty on the fit.

indicating that on both integrated and resolved scales, $\gamma$ is not closely tied to the distribution of old stars.

3. The behaviors of $U_{min}$ and $\bar{U}$ are very similar, as would be expected. In both of these parameters, we see that the slope of the resolved relationships is steeper than the equivalent integrated scaling relationship.

To summarize, the SED modeling results find a strong correlation between the minimum radiation field intensity, $U_{min}$, and the average radiation field intensity, $\bar{U}$, with quantities tied to SFR or sSFR. $U_{min}$ and $\bar{U}$ also show correlation on resolved scales with stellar mass surface density, $\Sigma_*$, suggesting that both older stars and recent star formation contribute to setting the typical radiation field intensity. The power-law component does not show significant correlation with $\Sigma_*$, but does with $\Sigma_{SFR}$, in agreement with the expectation that these more intense radiation fields are caused by proximity to star-forming regions.

The $\gamma$ parameter is weakly, positively correlated with $M_*$ for integrated galaxies, which is unexpected given that it is tied to higher radiation field intensities. The correlation we observe may be due to the fact that the power-law plus delta-function model applied to the entire galaxy enforces that a single $U_{min}$ must represent all areas of the galaxy. If $U_{min}$ is tracing $\Sigma_{SFR}$,

then the galaxy actually has a range of resolved $U_{min}$ values. The integrated $U_{min}$ may lie at the low end of the distribution of resolved $U_{min}$ values and push all higher $U_{min}$ values into the power-law distribution.

In general, $U_{min}$, $\bar{U}$, and $\gamma$ recovered by the power-law plus delta function model fits, particularly on resolved scales, appear in good agreement with expectations for the behavior of these components of the radiation field. This is discussed further in Section 6.1.

For predicting the resolved average radiation field intensity, $\bar{U}$, the best single quantity is $\Sigma_{SFR}$ (followed closely by $\Delta MS_R$). For $\gamma$, it is $\Delta MS_R$ and $\Sigma_{SFR}$. We note that $\Delta MS_R$, as defined in Equation (5), is not trivially applied outside our data set, so $\Sigma_{SFR}$ is likely the most accessible predictor for resolved radiation field quantities. In Section 5.6 we explore if the prediction of $\bar{U}$ is improved using a combination of $\Sigma_*$ and $\Sigma_{SFR}$.

### 5.3. Dust Properties' Scaling Relations: Integrated Fits

The results of integrated fits of the dust mass and PAH fraction are shown in Figure 6. Trends evident in $q_{PAH}$ in these plots include:
<figure><img src="image_11.png" /><figcaption>**Figure 7.** 2D histograms of $U_{\min}$ (top row), $\gamma$ (middle row), and $\bar{U}$ (bottom row) fit for individual lines of sight at the SPIRE 250 µm resolution, as a function of $sSFR_R$ (first column), $\Sigma_{SFR}$ (second column), $\Sigma_*$ (third column), and $\Delta MS_R$ (fourth column), color coded by the number of hits. The vertical white-dashed lines are the limits of completeness (Section 4.4). The white circles are the running medians, with error bars as one standard deviation, and the filled circles are those used to fit the scaling relations. We mark the slope in each panel.</figcaption></figure>

1. $q_{\text{PAH}}$ is not strongly correlated with SFR but is strongly, negatively correlated with sSFR and $\Delta$MS. This result suggests that at a fixed $M_*$, galaxies with higher SFR have lower PAH fractions. This may be tied to the destruction of PAHs in ionized gas, as discussed in Section 6.2.

2. $q_{\text{PAH}}$ is positively correlated with $M_*$. This scaling relation may result from a combination of a metallicity trend (lower $q_{\text{PAH}}$ at lower metallicity plus the mass-metallicity relationship) and an sSFR-related destruction of PAHs in H II regions (higher sSFR at lower $M_*$).

As discussed in Section 6.2, our results compiling a large sample of resolved nearby galaxies provide key context for previous studies of resolved PAH fraction in smaller samples. Our scalings also agree with higher-resolution PAH fraction mapping, which clearly shows the lower PAH fractions in H II regions likely due to destruction of PAHs (e.g., K. D. Gordon et al. 2008; D. Paradis et al. 2011; R. Paladini et al. 2012; M. Relaño et al. 2016; J. Chastenet et al. 2019, 2023; O. V. Egorov et al. 2023; J. Sutter et al. 2024). In Section 6.2, we discuss these results further and compare to past literature results on these correlations.

The bottom row of Figure 6 shows the integrated scaling relations of dust mass ($M_d$) as a function of SFR, $M_*$, sSFR, and $\Delta$MS. These panels show a few important trends:

1. $M_d$ is well correlated with both SFR and $M_*$. We find a remarkably linear slope in $M_d-M_*$ and small scatter. The correlation with SFR is slightly weaker and sublinear.

2. We find that $M_d$ is uncorrelated with sSFR, but shows a slight positive correlation with $\Delta$MS. This suggests that the lack of sSFR trend could result from competing effects of lower dust-to-gas ratios (DGRs) for lower $M_*$ galaxies and higher $M_{\text{gas}}/M_*$ for galaxies with positive offsets above the main sequence.

Trends of $M_d$ as a function of various galaxy-integrated quantities have been explored in detail in the literature. In Section 6.3, we discuss our results in comparison to past measurements and scaling relationships.

5.4. Dust Properties' Scaling Relations: Resolved Fits

In Figure 8 we show the resolved correlations between $\Sigma_d$ and $q_{\text{PAH}}$ as a function of $\Sigma_{SFR}$, $\Sigma_*$, $sSFR_R$, and $\Delta MS_R$. We note the following results about the behavior of $q_{\text{PAH}}$:

1. The resolved PAH fraction decreases steeply with $\Sigma_{SFR}$, $sSFR_R$, and $\Delta MS_R$. The steepest decrease occurs as a function of $sSFR_R$, where, over the span of slightly more than 1 dex in $sSFR_R$, the $q_{\text{PAH}}$ drops from values ~4%–5%,
<figure><img src="image_12.png" /><figcaption>Figure 8. Same as Figure 7, for fitted parameters $q_{\text{PAH}}$ (top row) and $\Sigma_d$ (bottom row).</figcaption></figure>

representative of the diffuse ISM of the Milky Way (B. T. Draine & A. Li 2007) to ~1% at log sSFR$_R$ ~−9. This trend is less steep in $\Delta\text{MS}_R$ and $\Sigma_{\text{SFR}}$, though both reach $q_{\text{PAH}}$ ~1% at their highest values. These trends suggest that the intensity of star formation in a given region is critical to setting the PAH fraction.

2. $q_{\text{PAH}}$ is weakly correlated with $\Sigma_*$, with a negative slope. This trend may result from the overall radial correlation of $\Sigma_*$ and $\Sigma_{\text{SFR}}$. The highest $\Sigma_*$ values tend to fall in galaxy centers, where $\Sigma_{\text{SFR}}$ can be high leading to lower $q_{\text{PAH}}$. The decrease at high $\Sigma_*$ could also be related to the less-efficient excitation of PAH vibrational modes where older stars contribute significantly to the radiation field (B. T. Draine et al. 2014, 2021; C. M. Whitcomb et al. 2024).

We note that this is the first compilation of resolved measurements of $q_{\text{PAH}}$ for a large, diverse sample of galaxies, where we can control for completeness effects in $\Sigma_{\text{SFR}}$ and $\Sigma_*$. Past studies of small samples have tended to select bright, highly star-forming galaxies that sit above the typical integrated or resolved star-forming main sequence. This makes interpreting trends as a function of other quantities like metallicity subject to biases due to incomplete sampling of sSFR. Our results clearly show a steep, negative dependence of $q_{\text{PAH}}$ on quantities related to SFR, which is likely related to their destruction in the ionized gas of H II regions. We discuss these insights in more detail in Section 6.2.

The bottom row of Figure 8 shows the resolved behavior of dust mass surface density, $\Sigma_d$. Here we note the following trends:

1. $\Sigma_d$ is positively correlated with $\Sigma_{\text{SFR}}$, $\Sigma_*$, and $\Delta\text{MS}_R$, but weakly, negatively correlated with sSFR$_R$.

2. The $\Sigma_d-\Sigma_{\text{SFR}}$ relationship can be interpreted as a reflection of the resolved Kennicutt–Schmidt relationship (F. Bigiel et al. 2008; A. K. Leroy et al. 2008), a topic we discuss further in Section 6.3.1.

3. $\Sigma_d$ is positively correlated with resolved main sequence offset, $\Delta\text{MS}_R$ and slightly negatively correlated with

$\text{sSFR}_R$. This reflects the fact that the resolved main sequence has a slope as a function of $\Sigma_*$. The trend of higher $\Sigma_d$ at a fixed $\Sigma_*$ for regions with higher $\Sigma_{\text{SFR}}$ (i.e., the positive correlation with $\Delta\text{MS}_R$) goes in the direction one would expect if those regions have higher gas surface densities. This agrees with resolved scaling relationships between gas, SFR, and stellar mass surface densities.

Comparing integrated and resolved results (Figure 9), the most dramatic differences are in the dependence of $q_{\text{PAH}}$ on $\Sigma_*$ and $\Sigma_{\text{SFR}}$. In the resolved measurements, $q_{\text{PAH}}$ decreases strongly with $\Sigma_{\text{SFR}}$, while the integrated $q_{\text{PAH}}$ is only weakly negatively correlated with SFR. We interpret this as a result of PAH destruction being related to the local intensity of star formation, which the resolved analysis better isolates. In the integrated analysis, $M_*$–metallicity and $M_*$–SFR scalings may overwhelm this signal.

### 5.5. Scaling Relations: Galaxy-averaged Surface Densities Fits

The results of the integrated analysis presented in Sections 5.1 and 5.3 are not straightforwardly comparable to the resolved measurements since they mix together “intensive” and “extensive” quantities, resulting in correlations that are more difficult to interpret. For example, $U_{\min}$ or $\bar{U}$ are intensive, average-like values, $q_{\text{PAH}}$ and $\gamma$ are fractions, while SFR and $M_*$ are extensive, sum-like integrated values. This means that a quantity like SFR can increase by adding together many low-SFR regions, while $\bar{U}$ or $q_{\text{PAH}}$ can stay the same under those conditions. This effect is likely to play a role in the different scalings we show in the resolved and integrated cases. To investigate this issue, we consider galaxy-averaged SFR, $\langle\Sigma_{\text{SFR}}\rangle$, and $M_*$, $\langle\Sigma_*\rangle$, as described in Section 3.8.2, to create intensive $M_*$ and SFR quantities by dividing by the effective area of the galaxy.

In Figure 10, we show the radiation field and dust parameters $U_{\min}$, $\bar{U}$, $\gamma$ and $q_{\text{PAH}}$ as a function of the galaxy-averaged SFR and $M_*$, noted $\langle\Sigma_{\text{SFR}}\rangle$ and $\langle\Sigma_*\rangle$. We also apply the same
<table><caption>**Table 4**
Coefficients to the Resolved Fits, $Y = aX + b$, Presented in Figures 7 and 8</caption><thead><tr><th rowspan="2">Y = Dust Parameter</th><th colspan="5">Binned Medians</th><th colspan="2">All Points</th></tr><tr><th>a ± ε<sub>a</sub><sup>fit</sup> ± ε<sub>a</sub><sup>bin</sup></th><th>b ± ε<sub>b</sub><sup>fit</sup> ± ε<sub>b</sub><sup>bin</sup></th><th>Pearson's ρ</th><th>RMSE</th><th>p-value</th><th>a'</th><th>p-value'</th></tr></thead><tbody><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(Σ<sub>SFR</sub>)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.42 ± 0.01 ± 0.01</td><td>1.08 ± 0.02 ± 0.02</td><td>0.99</td><td>0.235</td><td>&lt;0.001</td><td>0.4</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.39 ± 0.02 ± 0.01</td><td>-1.25 ± 0.04 ± 0.02</td><td>0.98</td><td>0.442</td><td>&lt;0.001</td><td>0.53</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.52 ± 0.02 ± 0.01</td><td>1.36 ± 0.03 ± 0.01</td><td>0.99</td><td>0.254</td><td>&lt;0.001</td><td>0.48</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-1.84 ± 0.13 ± 0.06</td><td>0.43 ± 0.21 ± 0.08</td><td>-0.96</td><td>1.498</td><td>&lt;0.001</td><td>-1.17</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Σ<sub>d</sub> [M<sub>☉</sub> pc<sup>-2</sup>])</td><td>0.49 ± 0.02 ± 0.01</td><td>-0.09 ± 0.03 ± 0.02</td><td>0.99</td><td>0.309</td><td>&lt;0.001</td><td>0.53</td><td>&lt;0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(Σ<sub>*</sub>)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.21 ± 0.01 ± 0.01</td><td>-0.21 ± 0.03 ± 0.03</td><td>0.98</td><td>0.264</td><td>&lt;0.001</td><td>0.22</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>-0.12 ± 0.01 ± 0.01</td><td>-1.82 ± 0.03 ± 0.03</td><td>-0.95</td><td>0.512</td><td>&lt;0.001</td><td>-0.08</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.21 ± 0.01 ± 0.01</td><td>-0.16 ± 0.03 ± 0.03</td><td>0.98</td><td>0.291</td><td>&lt;0.001</td><td>0.21</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-0.7 ± 0.09 ± 0.09</td><td>6.04 ± 0.26 ± 0.22</td><td>-0.92</td><td>1.657</td><td>&lt;0.001</td><td>-0.64</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Σ<sub>d</sub> [M<sub>☉</sub> pc<sup>-2</sup>])</td><td>0.29 ± 0.03 ± 0.04</td><td>-1.63 ± 0.07 ± 0.08</td><td>0.9</td><td>0.337</td><td>&lt;0.001</td><td>0.28</td><td>&lt;0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = log<sub>10</sub>(sSFR<sub>R</sub>)</td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.47 ± 0.06 ± 0.03</td><td>4.87 ± 0.59 ± 0.27</td><td>0.94</td><td>0.282</td><td>&lt;0.001</td><td>0.38</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.61 ± 0.06 ± 0.03</td><td>3.99 ± 0.56 ± 0.3</td><td>0.95</td><td>0.449</td><td>&lt;0.001</td><td>0.63</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.66 ± 0.09 ± 0.05</td><td>6.81 ± 0.89 ± 0.47</td><td>0.93</td><td>0.327</td><td>&lt;0.001</td><td>0.47</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-2.71 ± 0.33 ± 0.21</td><td>-22.9 ± 3.24 ± 2.13</td><td>-0.89</td><td>1.446</td><td>&lt;0.001</td><td>-2.57</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Σ<sub>d</sub> [M<sub>☉</sub> pc<sup>-2</sup>])</td><td>-0.15 ± 0.04 ± 0.06</td><td>-2.71 ± 0.35 ± 0.54</td><td>-0.72</td><td>0.377</td><td>0.01</td><td>-0.22</td><td>&lt;0.001</td></tr><tr><td colspan="8" style="text-align:center; font-weight:bold;">X = ΔMS<sub>R</sub></td></tr><tr><td>log<sub>10</sub>(U<sub>min</sub>)</td><td>0.54 ± 0.02 ± 0.01</td><td>0.12 ± 0.01 ± 0.01</td><td>0.99</td><td>0.239</td><td>&lt;0.001</td><td>0.5</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(γ)</td><td>0.61 ± 0.03 ± 0.02</td><td>-2.21 ± 0.02 ± 0.01</td><td>0.98</td><td>0.417</td><td>&lt;0.001</td><td>0.75</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Ū)</td><td>0.69 ± 0.03 ± 0.01</td><td>0.17 ± 0.02 ± 0.01</td><td>0.99</td><td>0.258</td><td>&lt;0.001</td><td>0.61</td><td>&lt;0.001</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>-2.31 ± 0.16 ± 0.13</td><td>4.41 ± 0.13 ± 0.06</td><td>-0.97</td><td>1.389</td><td>&lt;0.001</td><td>-1.94</td><td>&lt;0.001</td></tr><tr><td>log<sub>10</sub>(Σ<sub>d</sub> [M<sub>☉</sub> pc<sup>-2</sup>])</td><td>0.38 ± 0.02 ± 0.01</td><td>-1.21 ± 0.01 ± 0.0</td><td>0.98</td><td>0.347</td><td>&lt;0.001</td><td>0.43</td><td>&lt;0.001</td></tr></tbody></table>

**Note.** The quoted (a, b) values are for the displayed number of bins, which are running medians of the whole sample; the first uncertainty is the statistical error on these linear fits; the second uncertainty is the standard deviation on coefficients (a, b) when fitting binned data with a varying number of bins. Only the filled medians are considered for the fit, and chosen to have at least 250 points within the bin and above the completeness thresholds defined in the text (selection bias due to the 3σ dust emission fit). The Pearson's ρ coefficient is that of the fits shown. The RMSE is calculated on all dark-gray points. The last column shows a', the slope derived by fitting all points, without binning, and gives an additional sense of the uncertainty on the fit.

galaxy-averaged measurement to the total dust mass so that we look at the galaxy-averaged dust mass, $\langle \Sigma_d \rangle$. Here we see that most of the galaxy-averaged trends are very similar to the resolved trends, though the slopes are often slightly shallower. We report these slopes in Figure 9. Galaxy-averaged scalings may be of use particularly for distant galaxies, where estimates of $r_e$ are possible with HST and JWST imaging, but IR measurements may remain unresolved. Using the galaxy-averaged values provides a closer match to the resolved scaling behavior in most cases.

## 5.6. Predicting Radiation Field and Dust Properties

In the previous sections, we presented scaling relations between radiation field and dust parameters with single quantities related to stellar mass and SFR. For the purposes of predicting these properties, it is interesting to determine whether combinations of these variables yield improved predictions.

In Figure 11, we present 2D fits to a combination of $\Sigma_*$ and $\Sigma_{SFR}$ to determine a prediction for $\bar{U}$, $q_{PAH}$, and $\Sigma_d$. The results of these fits are also summarized in Table 5. In the top row of Figure 11, we show the 2D histograms between $\Sigma_{SFR}$ and $\Sigma_*$,

color-coded by the average $\bar{U}$, $q_{PAH}$, and $\Sigma_d$ value of all of the points that fall in that bin. The bins that are outside of the completeness bounds are shown with increased transparency and are not included in our fits. We fit a 2D plane of the form $Y = a \log_{10}(\Sigma_*) + b \log_{10}(\Sigma_{SFR}) + c$, and the results of these fits are listed in Table 5. In all cases, using two variables improves the resulting scatter compared to the single-parameter fits listed in Table 4.

In the left panels, we see that the prediction of $\bar{U}$ using both $\Sigma_{SFR}$ and $\Sigma_*$ results in RMSE of 0.2 dex, improving by 0.05 dex from the best single-parameter predictions (with $\Sigma_{SFR}$ alone). The diagonal color gradient shows that both $\Sigma_{SFR}$ and $\Sigma_*$ are correlated with $\bar{U}$.

In the center column of Figure 11, the top row shows a similar diagonal color gradient, implying that both $\Sigma_*$ and $\Sigma_{SFR}$ also affect the PAH fraction; however, given the coefficients of the 2D fit, $\Sigma_{SFR}$ matters more than $\Sigma_*$. Additionally, in this case, the contribution of $\Sigma_*$ appears to be positive, as opposed to being slightly negative in the 1D fit. This is likely the result of the $\Sigma_*-\Sigma_{SFR}$ information allowing separation of radial trends in these quantities. In general, the scatter in the predicted $q_{PAH}$ values is only improved a small amount by the 2D fits. For instance, as listed in Table 4, $\Delta\text{MS}_R$ alone can predict $q_{PAH}$
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Comparison of the slopes (a) of resolved (black symbols) and integrated (gray symbols) fits presented in Figures 5–8 and Tables 3 and 4. We also show the fit to $\langle$SFR$\rangle$ and $\langle M_*$$, galaxy-averaged values, shown in Figure 10. Here, each fit follows the form $Y = aX + b$, where we show only a in this figure. The open symbols mark where the high p-value signals no significant correlation between the axes. The behavior of the slopes of the $U_{\min}$ and $\bar{U}$ relationships are very similar, since under most circumstances the radiation field is dominated by the minimum radiation field.</figcaption></figure>

to RMSE of 1.39%. The best two-parameter fit ($\Sigma_{*}$ and $\Sigma_{SFR}$) improves this slightly to an RMSE of 1.2%.

We perform the same 2D fits to predict the dust surface density. In the one-parameter fits, $\Sigma_{SFR}$ provided the best RMSE of 0.31 dex. The two-parameter fits improve this slightly to 0.26 dex for the combination of $\Sigma_{*}$ and $\Sigma_{SFR}$. It is worth noting that these fits show that it is possible to predict the $\Sigma_{d}$ to within less than a factor of 2 scatter using only information about $\Sigma_{*}$ and $\Sigma_{SFR}$.

# 6. Discussion: Dust Properties across the Local Universe

## 6.1. Drivers of the Radiation Field Distribution in Galaxies

The radiation field intensity distribution in galaxies is of fundamental importance to understanding their dust emission. Radiation field intensity also plays critical roles in the physics of the ISM (regulating phase balance, H$_{2}$ formation, etc.). At the scales we resolve in nearby galaxies (hundreds of parsecs to kiloparsecs), there is inevitably a distribution of radiation field intensities heating the dust. Knowledge of this radiation field

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Galaxy-averaged SFR and $M_*$, $\langle$SFR$\rangle$ (left) and $\langle M_*$$\rangle$ (right), variations with the dust parameters $U_{\min}$, $\gamma$, $\bar{U}$, $q_{PAH}$, and galaxy-averaged dust mass, $\langle M_d \rangle$. The symbols show the binned medians and lines are linear fits to the filled medians, similarly to Figures 5–8.</figcaption></figure>

distribution is critical for extracting information about changes in dust properties from the SED. For instance, single-temperature modified blackbody models that fit an emissivity slope $\beta$ are subject to biases if the peak of the thermal SED is broadened by a radiation field distribution (L. K. Hunt et al. 2015). In addition, measuring changes in “very small grain” abundance is often degenerate with the assumed radiation field distribution (F. Galliano 2022).

Inferring the properties of this radiation field distribution can be approached in several ways. At the most detailed level, radiative transfer modeling of galaxies can attempt a full solution at some spatial resolution (e.g., I. De Looze et al. 2014; A. Nersesian et al. 2020). However, radiative transfer modeling with sufficiently high resolution and realistic dust, SFR, and $M_*$ distributions remains computationally challenging. Radiative transfer models for “pieces” of galaxies (K.-H. Law et al. 2018) are an alternative that yield insights
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Two-dimensional fits of $\bar{U}$, $q_{PAH}$, and $\Sigma_d$ as a function of a combination of the resolved stellar mass and star formation rate. Top row: 2D histograms of the $\Sigma_\star-\Sigma_{SFR}$ plane, where the bin values show the average $\bar{U}$ (left column), $q_{PAH}$ (center column), and $\Sigma_d$ (right column). The light bins are calculated for all pixels passing the $3\sigma$ S/N cut, and the bright bins are calculated for pixels passing that same cut and that have $\Sigma_\star$ and $\Sigma_{SFR}$ above the completeness thresholds. The 2D fit uses the standard deviations within each bin as errors in the linear regression. Bottom row: the x-axis shows the best-fit using both $\Sigma_\star$ and $\Sigma_{SFR}$. The y-axis shows the observed data. The best-fit coefficients are reported in Table 5.</figcaption></figure>

<table><caption>**Table 5**
Coefficients of Plane Fits: $Y = a \times \log_{10}(\Sigma_\star) + b \times \log_{10}(\Sigma_{SFR}) + c$</caption><thead><tr><th>Y = Dust Parameter</th><th>a</th><th>b</th><th>c</th><th>RMSE</th></tr></thead><tbody><tr><td>log<sub>10</sub>(Ū)</td><td>-0.22 ± 0.02</td><td>0.67 ± 0.02</td><td>2.14 ± 0.06</td><td>0.200</td></tr><tr><td>q<sub>PAH</sub> [%]</td><td>1.25 ± 0.09</td><td>-2.44 ± 0.09</td><td>-3.68 ± 0.29</td><td>1.196</td></tr><tr><td>log<sub>10</sub>(Σ<sub>d</sub> [M<sub>☉</sub> pc<sup>-2</sup>])</td><td>0.41 ± 0.03</td><td>0.22 ± 0.02</td><td>-1.42 ± 0.08</td><td>0.257</td></tr></tbody></table>

into the radiation field distribution. A less computationally intensive approach that provides information on radiation field distributions is “energy-balance” SED fitting, using codes like CIGALE (S. Noll et al. 2009; M. Boquien et al. 2019), Prospector-$\alpha$ (J. Leja et al. 2017), or MAGPHYS (E. da Cunha et al. 2010). These codes use assumptions about how attenuated starlight from various stellar populations translates to the reradiated IR SED. Fitting the UV to IR SED with these models provides constraints on the radiation field distributions, given those assumptions. Such modeling has been widely used for unresolved galaxies (e.g., L. K. Hunt et al. 2019; A. Nersesian et al. 2019) and has been applied to small sets of resolved galaxies as well (e.g., M. Boquien et al. 2012; M. Decleir et al. 2019; Abdurro’uf et al. 2022), though the spatial scales on which dust energy balance assumptions break down have not been fully quantified.

Another approach, and the one we use in this work, is to assume an empirically motivated shape for the radiation field distribution and constrain its parameters based on comparison to observed mid- to far-IR SEDs. This is possible because, having fixed the grain population properties, the radiation field intensity sets the equilibrium dust temperature and determines

the intensity of the mid-IR stochastic dust emission (which is ~linear with $U$). A distribution of radiation field intensities will therefore broaden the thermal peak of the SED while maintaining the linear proportionality of the dust-mass weighted $\bar{U}$ and the stochastic emission. We have adopted a simple model to describe the radiation field distribution, following B. T. Draine & A. Li (2007), with a power-law plus delta function. The delta function is intended to describe dust in the diffuse ISM heated by the average interstellar radiation field. The power-law distribution is intended to describe the higher-intensity radiation fields in the vicinity of star-forming regions. While we do not have information about the true distribution of radiation fields in our galaxies to compare against, we can compare the resulting correlations of the radiation field parameters with environment to assess their plausibility.

In Sections 5.1 and 5.2, we described the integrated and resolved correlations of radiation field distribution with galaxy properties. We found that overall, the average radiation field intensity, $\bar{U}$, was similar to $U_{min}$, and the fraction of the dust luminosity from the power-law distribution, $\gamma$, was small. This suggests that the majority of dust heating is due to the average
interstellar radiation field, not to the very high-intensity regions
near current star formation. We also investigated how $\bar{U}$ was
correlated with galaxy properties. From our resolved analysis,
we found that $\bar{U}$ depends on both $\Sigma_{SFR}$ and $\Sigma_*$, suggesting both
current star formation and older stellar populations contribute
to the average radiation field intensity. These results are in good
agreement with radiative transfer studies, which suggest that
older stellar populations play a significant role (~ tens of %;
I. De Looze et al. 2014; S. Viaene et al. 2017; A. Nersesian et al.
2020; K.-H. Law et al. 2021; M. T. Rushton et al. 2022) in dust
heating. This is also the case for results of UV-IR energy
balance SED modeling (e.g., CIGALE), which find important
contributions from older stellar populations to dust heating
(A. Nersesian et al. 2019).

One clear trend from our analysis is the lack of correlation of $\gamma$ (the fraction of the dust luminosity from the power-law component) with quantities describing the stellar mass ($M_*$ or $\Sigma_*$) and a strong correlation with parameters describing star formation ($\Sigma_{SFR}$, sSFR$_R$, $\Delta MS_R$). This can be seen in Figures 5 and 7. This behavior for the power-law component is expected, and suggests the SED modeling is doing a reasonable job capturing the high-intensity parts of the radiation field.

In total, the average radiation field, $\bar{U}$, in a ~kiloparsec region of a galaxy has contributions from recent star formation both through $\gamma$ and through $U_{\min}$. In Section 5.6, we showed that plane fits to a combination of $\Sigma_{SFR}$ and $\Sigma_*$ allow $\bar{U}$ to be predicted to ~0.2 dex accuracy. This behavior suggests that dust heating from young stellar populations is strongly correlated with sSFR. These observations agree with the results of radiative transfer models, which also find a strong correlation between the young stellar heating fraction and sSFR (e.g., A. Nersesian et al. 2020), and with the results of energy balance SED fitting (M. Boquien & S. Salim 2021), which suggest a stronger trend between the dust temperature (here equivalent to $\bar{U}$) and the sSFR, rather than SFR or $M_*$ alone. A strong sSFR dependence is also consistent with analytical modeling results from H. Hirashita & I.-D. Chiang (2022) and I.-D. Chiang et al. (2023). The correlation of dust temperature (and $\bar{U}$) with sSFR on galaxy-integrated scales is a topic of great interest in studies of high-redshift galaxies, a topic we will revisit in Section 6.5.

## 6.2. The Dependence of PAH Fraction on Integrated and Resolved Galaxy Properties

PAHs are a critical component of the ISM, participating in
many processes like photoelectric heating, H₂ formation, and
UV attenuation (A. G. G. M. Tielens 2008). However,
fundamental aspects of their life cycle remain mysterious.
Two key observations that provide insights into the PAH life
cycle include: (1) their abundance relative to dust mass strongly
depends on metallicity (C. W. Engelbracht et al. 2005, 2008;
S. C. Madden et al. 2006; B. T. Draine et al. 2007;
K. D. Gordon et al. 2008; K. M. Sandstrom et al. 2010;
A. Rémy-Ruyer et al. 2015; G. Aniano et al. 2020; I. Shivaei
et al. 2024; C. M. Whitcomb et al. 2024), and (2) they are
destroyed in H II regions (D. Cesarsky et al. 1996; M. S. Povich
et al. 2007; J. Chastenet et al. 2019, 2023; O. V. Egorov et
2023; J. Sutter et al. 2024).

Observations of the Magellanic Clouds provide particularly
clear evidence for both of these trends. J. Chastenet et al. (2019)
found that at the LMC's metallicity of 1/2 $Z_\odot$, the diffuse neutral
gas had PAH fractions similar to the Milky Way diffuse neutral

gas (~4%–5%; B. T. Draine & A. Li 2007), while H II regions
showed steep drops in PAH fraction. Moving to the SMC, at
1/5 $Z_\odot$, the PAH fraction in the diffuse ISM was far lower
(~1%), with H II regions showing yet lower values. These
observations suggest a steep metallicity dependence in PAH
fraction that sets in below ~1/2 $Z_\odot$. Recent work by
C. M. Whitcomb et al. (2024) measured $q_{PAH}$ as a function of
radius in M101, NGC 628, and NGC 2403, revealing a steep
drop in at a threshold metallicity of 2/3 $Z_\odot$. With the mid-IR
resolution of JWST, it is now possible to separate H II regions
from the diffuse ISM in many nearby galaxies and examine
trends in mid-IR traced PAH fraction both in and out of
H II regions. J. Sutter et al. (2024) presented such observations
for 19 nearby galaxies, showing an approximately constant PAH
fraction outside H II regions and steep drops within them at
~ $Z_\odot$. They demonstrate that ≳kiloparsec-scale PAH fraction
in ~ $Z_\odot$ galaxies is highly correlated with sSFR, since it traces a
relative fraction of the ISM that is currently in H II regions.

Our observations provide a comprehensive view of the
~kiloparsec scale and integrated $q_{PAH}$ in nearby galaxies and
allow us to control for completeness in measuring the scaling of
$q_{PAH}$ with environment. One of our key findings is that $q_{PAH}$ is a
steeply decreasing function of sSFR and sSFR$_R$ (Figures 6 and
8): $q_{PAH}$ is strongly, negatively correlated on resolved scales
with all quantities that trace the local intensity of star formation
($\Sigma_{SFR}$, sSFR$_R$, and $\Delta MS_R$). On galaxy-integrated scales, strong
negative correlations with sSFR and $\Delta MS$ persist, though the
$q_{PAH}$ does not vary strongly with galaxy-averaged SFR. It is
clear from this summary that the local intensity of star formation
is a primary agent in setting the PAH fraction in galaxies, likely
due to the destruction of PAHs in H II regions.

Strong trends in PAH fraction correlated with sSFR have also
been seen in integrated galaxy measurements (A. Rémy-Ruyer
et al. 2015; A. Nersesian et al. 2019; F. Galliano et al. 2021) and
resolved studies (e.g., the trends with 70 µm surface brightness in
G. Aniano et al. 2020). Using CIGALE (M. Boquien et al. 2019),
and The Heterogeneous dust Evolution Model for Interstellar
Solids (THEMIS; A. P. Jones et al. 2017), A. Nersesian et al.
(2019) fit integrated UV-to-far-IR SEDs of the DustPedia sample
(J. I. Davies et al. 2017; C. J. R. Clark et al. 2018). They studied
the variations of the small carbonaceous grain content as a
function of galaxy morphology and found that the fraction
of small grains decreases with sSFR after a threshold
$\log_{10}(\text{sSFR}) \ge -10.5$, in good agreement with our findings.
They also suggest that this decrease of the small grain fraction is
due to the harder radiation fields found at higher sSFR. A. Rém-
y-Ruyer et al. (2015) found similar trends using the KINGFISH
and Dwarf Galaxy Survey (S. C. Madden et al. 2013) galaxy
samples. Our findings align very well with the previously
observed scaling relations, but now extend them to resolved
scales across the galaxy population.

Because of the shape of the star-forming main sequence,
sSFR is higher for lower $M_*$ galaxies (and sSFR$_R$ is higher for
lower $\Sigma_*$). To attempt to separate out $M_*$ trends (which may be
tied to metallicity via the mass–metallicity relationship;
C. A. Tremonti et al. 2004), we investigated the correlation
of PAH fraction with offset from the main sequence ($\Delta MS$ and
$\Delta MS_R$). The observed trends of $q_{PAH}$ with $\Delta MS$ and $\Delta MS_R$
are shallower than those with sSFR and sSFR$_R$, suggesting that
sSFR trends may include a combination of both metallicity and
sSFR effects. However, we also observe a slightly negative
slope of $q_{PAH}$ as a function of $\Sigma_*$. The observed dependence of
$q_{\text{PAH}}$ on metallicity, plus the generally increasing metallicity toward the centers of galaxies would lead us to expect a positive trend between $q_{\text{PAH}}$ and $\Sigma_{*}$. Our observed trend, which we note does have significant scatter, may be the result of the higher $\Sigma_{\text{SFR}}$ and sSFR$_R$ toward galaxy centers dominating over the gradient of metallicity and $\Sigma_{*}$ in setting $q_{\text{PAH}}$. Another option may be the contribution of softer radiation fields to heating PAHs at high $\Sigma_{*}$. Such an effect has been found to alter $q_{\text{PAH}}$ derived from mid- to far-IR modeling in the center of M31 (B. T. Draine et al. 2014) and in other galaxy centers (C. M. Whitcomb et al. 2024).

It is worth noting that in many studies of the metallicity dependence of the PAH fraction, galaxy samples have inevitably been limited by the need to find bright, low-metallicity targets. Selecting low-metallicity galaxies with detectable mid-IR emission has favored galaxies with unusually high SFRs for their masses (e.g., “blue compact dwarfs” like I Zw 18, Henize 2-10, Haro 11, and others). Given the strong dependence of PAH fraction on sSFR and $\Delta$MS, focusing on these highly star-forming galaxies may lead to overestimation of metallicity effects. This is also true for targeted spectroscopic studies of individual star-forming regions in galaxies—without sampling “average” star formation intensities, these observations may lead to biased conclusions about the strength of the metallicity dependence of the PAH fraction. We do not have detailed metallicity information for most galaxies in our sample, and our sample is not complete down to very low metallicity (stellar mass), so we cannot draw strong conclusions about metallicity dependence for $q_{\text{PAH}}$. Our strong trends of $q_{\text{PAH}}$ with $\Delta$MS and $\Delta$MS$_R$ show that there are large variations in $q_{\text{PAH}}$ even controlling for metallicity. Separating sSFR and metallicity effects therefore requires care toward completeness and unbiased samples of galaxies. Future studies with JWST, given its high sensitivity, should be able to disentangle the dependence of PAH fraction on sSFR and metallicity, given careful sample selection.

### 6.3. The Relationship between Dust Mass, Stellar Mass, and Star Formation Rate

Trends in $M_d$ with SFR and $M_*$ likely reflect the combined effects of several well-known scaling relationships: (1) the mass–metallicity relationship or fundamental mass–metallicity–SFR relationship (R. Maiolino & F. Mannucci 2019), which shows that lower $M_*$ galaxies tend to have lower metallicity; (2) the Kennicutt–Schmidt relationship (R. C. Kennicutt & N. J. Evans 2012), which describes the scaling between gas mass and SFR; (3) the gas mass to stellar mass scaling (A. Saintonge & B. Catinella 2022), which describes the behavior of gas fraction over the star-forming main sequence; and (4) the dependence of the dust-to-metals ratio on metallicity (which sets the mapping between metallicity and DGR; A. Rémy-Ruyer et al. 2014; P. De Vis et al. 2019; C. Péroux & J. C. Howk 2020). From the combination of these scaling relationships, we expect predictable trends in $M_d$ to arise over the SFR–$M_*$ space. Such trends have been explored by a number of studies for integrated measurements of galaxies (L. Cortese et al. 2012; P. De Vis et al. 2017a; V. Casasola et al. 2020; I. De Looze et al. 2020). These scaling relationships also have resolved behavior (F. Bigiel et al. 2008; A. K. Leroy et al. 2008; S. F. Sánchez et al. 2014; B. C. Hsieh et al. 2017; S. L. Ellison et al. 2018; L. Lin et al. 2019), which should similarly result in correlations of $\Sigma_d$ with $\Sigma_{\text{SFR}}$, $\Sigma_*$, sSFR$_R$, and $\Delta$MS$_R$. Our observations provide a powerful database with which to explore

these trends, carefully controlling for completeness in $\Sigma_d$ and selecting where we should be able to measure accurate SFRs.

On integrated scales, above our sSFR cut ($\log_{10}(\text{sSFR}) < -11$), we find a remarkably constant $M_d/M_*$ ratio spanning $\log_{10}(M_*) = 9-11$ (indicated by the slope of ~1 in Figure 6 and the nearly constant ratio in the bottom-left panel of Figure 12). A variety of other integrated studies of galaxies in the local Universe have drawn different conclusions, finding negative trends of $M_d/M_*$ with increasing $M_*(L. Cortese et al. 2012; P. De Vis et al. 2017a; V. Casasola et al. 2020; I. De Looze et al. 2020). I. De Looze et al. (2020) found a slope of -0.22 fitting the integrated SEDs of 423 galaxies, a combination of the JINGLE (A. Saintonge et al. 2018), KINGFISH (R. C. Kennicutt et al. 2011), HAPLESS (C. J. R. Clark et al. 2015), HiGH (P. De Vis et al. 2017a), and some of the HRS (A. Boselli et al. 2010) samples. We can reproduce negative slopes by relaxing our sSFR cut, as can be seen in Figure 12. Here, the light-gray points show galaxies that were excluded by our sSFR cut—they systematically fall below the average $M_d/M_*$. Including these points in our fit would yield a negative slope for $M_d/M_*$ versus $M_*$. The primary driver of this trend is likely to be the decreased gas-to-$M_*$ fraction in these low sSFR galaxies (L. Cortese et al. 2012; A. Saintonge & B. Catinella 2022).

Due to the mass–metallicity relationship, lower $M_*$ galaxies are expected to have lower metallicity, and hence lower DGR. In addition, lower-mass galaxies are observed to have higher gas fractions and higher sSFR. These DGR and gas fraction trends should therefore compete with each other in setting the $M_d/M_*$ scaling as a function of $M_*$. In Figure 12, we show in the bottom-right panel, the sSFR–$M_*$ space color-coded by $M_d/M_*$. The figure shows that on or near the star-forming main sequence, the trends at least approximately cancel out, leaving a ~ fixed $M_d/M_*$ ratio >2 orders of magnitude in $M_*$.

Distinct from the constant $M_d/M_*$ ratio, we see for integrated galaxies, on resolved scales $\Sigma_d/\Sigma_*$ is a decreasing function of $\Sigma_*$, with slope -0.51, as shown in Figure 12 in the top-left panel. This trend is very similar to the slope of -0.56 observed in a sample of 18 face-on spiral galaxies from DustPedia (V. Casasola et al. 2022) and also resembles high spatial resolution measurements of the $\Sigma_d-\Sigma_*$ scaling in M31 (S. Viaene et al. 2014). This measurement indicates that there is a sublinear power-law slope for the $\Sigma_d-\Sigma_*$ relationship (also shown in Figure 8). Furthermore, because the majority of galaxies have negative metallicity gradients, a given $\Sigma_d$ value drawn from the inner or central part of a galaxy will tend to correspond to less $\Sigma_{\text{gas}}$ than the same $\Sigma_d$ in the outskirts of the galaxy. This implies that most galaxies will also have a sublinear $\Sigma_{\text{gas}}-\Sigma_*$ relationship.

In outer disks, where H I represents most of the ISM mass, such a sublinear $\Sigma_{\text{gas}}-\Sigma_*$ relationship is in good agreement with the frequent observation that the atomic gas disks of galaxies often extend to larger radii compared to the stellar disks, with the $\Sigma_{\text{gas}}/\Sigma_*$ ratio increasing with increasing galactocentric radius (e.g., F. Bigiel et al. 2010; R. C. Kennicutt & N. J. Evans 2012; J. Wang et al. 2016). However, at $\Sigma_* \gtrsim 100 M_\odot \text{ pc}^{-2}$, where the ISM is H$_2$-dominated (A. K. Leroy et al. 2008; A. Schruba et al. 2011), many observations show a ~linear power-law relationship between CO emission and $\Sigma_*(A. K. Leroy et al. 2008; A. Schruba et al. 2011; L. Lin et al. 2019; S. F. Sánchez 2020; I. Pessa et al. 2022), suggesting a linear relationship between $\Sigma_{\text{mol}}-\Sigma_*$. This difference in the CO-$\Sigma_*$ and $\Sigma_d-\Sigma_*$ slopes may be the result of systematic variations in the CO-to-H$_2$ conversion factor ($\alpha_{\text{CO}}$). Recent work
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Left column: specific dust mass as a function of stellar mass on resolved (top) and integrated (bottom) scales. The light-gray points are those with $log_{10}(sSFR) < -11$. The medians are measured on the dark-gray points. Error bars show the scatter within the bin; the solid black line is the best-fit using the filled medians, whose slope is reported in the label; the first uncertainty is the statistical error on the fits, and the second uncertainty is the standard deviation on coefficients when fitting binned data with a varying number of bins. Orange dashed lines show the best fit from the literature: the best-fit from V. Casasola et al. (2022) for a few galaxies in the DustPedia sample, and from I. De Looze et al. (2020) for a large sample of integrated measurements. Right column: main sequence color-coded by the specific dust mass. In the resolved panel (top), the transparency shows all pixels, while the opaque bins cut out pixels with $log_{10}(sSFR) < -11$ and above the completeness threshold in $\Sigma_{*}$. The line shows our main-sequence fit from Section 3.8.3. The bottom panel shows the integrated values from this work for the main sequence, and the white-dashed line shows the main-sequence fit by A. K. Leroy et al. (2019).</figcaption></figure>

by I.-D. Chiang et al. (2024) used measurements of CO, H I, $\Sigma_{d}$, and metallicity for a sample of 37 nearby galaxies, with an assumption of a constant dust-to-metals ratio, to infer $\alpha_{CO}$ that decreases with increasing $\Sigma_{*}$, $\alpha_{CO}^{2-1} \propto \Sigma_{*}^{-0.5}$ (for CO $J=2-1$, more shallowly for CO $1-0$). From our compilation of 819 resolved galaxies, we can see that the shallow $\Sigma_{d}-\Sigma_{*}$ is a general feature of the sample. This observation, combined with the widely observed, linear CO-$\Sigma_{*}$ relationship, suggests that a decreasing $\alpha_{CO}$ with $\Sigma_{*}$ is a general trend across the galaxy population.

As in the integrated case, we find that $\Sigma_{d}/\Sigma_{*}$ is positively correlated with the offset above the resolved main sequence ($\Delta MS_{R}$), as shown in Figure 12 in the top-right panel. Assuming that fixed $\Sigma_{*}$ controls for metallicity and DGR

variations to first order, this trend would suggest that points are offset above the resolved main sequence (high $\Delta MS_{R}$) because they have a larger dust (and gas) surface density. This result agrees with recent work showing that the correlation properties of $\Sigma_{mol}-\Sigma_{SFR}-\Sigma_{*}$ suggests that offset from the resolved main sequence is primarily driven by having higher $\Sigma_{mol}$ (W. M. Baker et al. 2022). We note, however, that conclusions drawn about offsets in $\Sigma_{mol}$ are subject to uncertainties in $\alpha_{CO}$, while $\Sigma_{d}$ may have a less variable relationship to gas.

### 6.3.1. The $\Sigma_{d}-\Sigma_{SFR}$ Relationship

The $\Sigma_{d}-\Sigma_{SFR}$ relationship is closely tied to the resolved star formation scaling law or Kennicutt–Schmidt relationship, which
relates gas surface density ($\Sigma_{gas}$) and $\Sigma_{SFR}$ (R. C. Kennicutt & N. J. Evans 2012), as well as the closely related concept of the local gas depletion time, $\tau_{dep} = \Sigma_{gas}/\Sigma_{SFR}$. There has been a substantial effort to characterizing this relationship. The slope is found to vary depending on whether the ISM is H I- or H$_2$-dominated (F. Bigiel et al. 2008; A. K. Leroy et al. 2008). Shorter gas depletion times and more nearly linear slopes are often found in regions of galaxies dominated by molecular gas, while the H I-dominated parts of galaxies show a much wider range of depletion times, often including very long $\tau_{dep}$ in outer galaxies, and steeper $\Sigma_{SFR}-\Sigma_{gas}$ relationships. This is understood to reflect that stars form from molecular gas, and that the balance of phases in the ISM is tied to local environmental conditions, including dynamical equilibrium pressure in the ISM and the DGR (L. Blitz & E. Rosolowsky 2006; A. K. Leroy et al. 2008; J. Sun et al. 2020; C. Eibensteiner et al. 2024).

The measured slope in the H$_2$-dominated ISM is sensitive to the assumed behavior of the CO-to-H$_2$ conversion factor (J. Sun et al. 2023; Y.-H. Teng et al. 2024). However, many measurements of the CO-to-H$_2$ conversion factor ($\alpha_{CO}$) are based on $\Sigma_d$, with built-in assumptions about the DGR. This makes all aspects of this problem closely connected (I.-D. Chiang et al. 2024). For high-redshift galaxies, the $\Sigma_d-\Sigma_{SFR}$ scaling may be more observationally accessible than the relationship between H I or H$_2$ and $\Sigma_{SFR}$, making our observed relation useful in its own right. Since the ratio $\Sigma_d/\Sigma_{SFR}$ may be even more accessible at high-z than the scaling, we note that for an intermediate $\Sigma_d = 0.1 M_{\odot} pc^{-2}$, our fit implies a “dust depletion time” of

$$ \tau_{\text{dust dep}} = \frac{\Sigma_{\text{dust}}}{\Sigma_{\text{SFR}}} \approx 8.7 \text{ Myr}. \qquad (8) $$

Though star formation is clearly fueled by gas, not only dust, this ratio expressed the SFR per unit dust and may be an interesting point of comparison for high-z systems or other populations.

As has been extensively discussed in the literature (G. A. Blanc et al. 2009; A. K. Leroy et al. 2013; R. Shetty et al. 2014), methodological choices in binning and power-law fitting play an important role in determining the resulting slope of the star formation scaling relation. In Figure 13, we show slopes for fits to the binned averages, with binning set either in x or y, in light orange. Because of the intrinsic scatter in the data, binning by one or the other axis changes the slope and biases the derived scaling relation. To minimize these biases related to the choice of the dependent variable, we use the “bisector” fit (in black), which takes the mean of the two slopes. The fit is described by the equation:

$$ \log_{10}(\Sigma_{\text{SFR}}) = (1.40 \pm 0.08) \times \log_{10}(\Sigma_d) - (0.54 \pm 0.10), \qquad (9) $$

with $\Sigma_{SFR}$ in $M_{\odot} kpc^{-2} yr^{-1}$ and $\Sigma_d$ in $M_{\odot} pc^{-2}$. The slope we fit is in good agreement with recent work by I.-D. Chiang et al. (2023), who also derived a $\Sigma_{SFR}-\Sigma_d$ relationship.

Our observed $\Sigma_d-\Sigma_{SFR}$ fit, assuming a constant DGR, provides a measurement of the $\Sigma_{SFR}-\Sigma_{gas}$ scaling law. Our observed slope of 1.4 is slightly steeper than the slope of 1.2 obtained for only molecular by J. Sun et al. (2023) using dust-based $\alpha_{CO}$ estimates (their B13 values) or similar results by Y.-H. Teng et al. (2024). This likely reflects that our measurements include dust and SF associated with both the inner H$_2$-dominated and the outer H I-dominated parts of disks. At DGR = 0.01, typical for solar metallicity ISM, the lowest surface densities included in our fit,

<figure><img src="image_17.png" /><figcaption>Figure 13. Two-dimensional histogram and scaling relation of the dust surface density, $\Sigma_d$, and the star formation rate surface density, $\Sigma_{SFR}$. The orange medians and lines are derived from binning either in x or y, and the black line shows the mean of these slopes, with the reported equation.</figcaption></figure>

$\log_{10}(\Sigma_d [M_{\odot} pc^{-2}]) \sim -1.5$, would correspond to a total gas mass surface density as low as $\Sigma_{gas} \sim 3 M_{\odot} pc^{-2}$. This gas surface density range extends well into the H I-dominated parts of disks, showing that we are not only sampling the more linear $\Sigma_{SFR}-\Sigma_{gas}$ relationship in H$_2$-dominated regions. This appears to be in good agreement with measurements of the $\Sigma_{SFR}-\Sigma_{gas}$ scaling relation that include both H I and H$_2$ (F. Bigiel et al. 2008; A. Schruba et al. 2011; R. C. Kennicutt & N. J. Evans 2012).

Alternatively to fitting the $\Sigma_d-\Sigma_{SFR}$ slope, we could assume a gas depletion time to convert the observed $\Sigma_{SFR}$ into a gas mass surface density, from which we could determine the DGR (i.e., “inverting the Kennicutt-Schmidt relation”; e.g., R. Genzel et al. 2012; A. Schruba et al. 2012). One option would be to estimate a $\tau_{dep}$ that depends on $\Sigma_{gas}$ according to the galaxy-integrated $\Sigma_{gas}-\Sigma_{SFR}$ relation from R. C. Kennicutt (1998), which has a power-law slope of 1.4 and intercept of $2.5 \times 10^{-4} M_{\odot} yr^{-1} kpc^{-2}$. In this case, converting $\Sigma_{SFR}$ to $\Sigma_{gas}$ and comparing to $\Sigma_{dust}$, we estimates a median DGR of $\log_{10}(\text{DGR}_{KS}) = -2.2 \pm 0.3$ across our whole data set. This value is slightly lower, but relatively consistent with the range of $M_d/M_H$ ratios found in Milky Way (i.e., solar metallicity) dust models: 0.01 in the M. Compiègne et al. (2011) model, 0.0074 in THEMIS (A. P. Jones et al. 2017), or 0.01 in B. S. Hensley & B. T. Draine (2023). Additionally, we make the same measurement with the $\Sigma_{H_2}-\Sigma_{SFR}$ relationship from J. Sun et al. (2023), appropriate for the H$_2$-dominated ISM. We use their “mKS” prescription and find $\log_{10}(\text{DGR}_{mKS}) = -2.4 \pm 0.4$. We note that J. Sun et al. (2023) adopted a metallicity-varying CO-to-H$_2$ conversion factor, which is calibrated using dust as a tracer of total gas (A. D. Bolatto et al. 2013). In fact, it is worth noting that observations of DGR, the CO-to-H$_2$ conversion factor, and molecular gas depletion time ($\Sigma_{mol}/\Sigma_{SFR}$) are all observationally connected through the uncertainty in converting CO to H$_2$.

## 6.4. $q_{PAH}$ and WISE Bands

Several studies have focused on using band ratios of mid-IR photometry as a tracer of PAH fraction, often contrasting a PAH-dominated band with a continuum-dominated one. This was done with Spitzer at ~8 and 24 µm (e.g., C. W. Engelbracht et al. 2005, 2008; J. D. T. Smith et al. 2007; A. R. Marble et al. 2010; K. V. Croxall et al. 2012), and more recently with JWST (e.g., J. Chastenet et al. 2023; O. V. Egorov et al. 2023;
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Fitted $q_{PAH}$ from the resolved analysis as a function of $log_{10}$(WISE 3/WISE 4) ratio. We fit a line through the bins that pass the same cuts used in other figures, for consistency, marked in black-filled symbols. For prediction purposes only, we show an additional fit in the form of a quadratic fit (see the text).</figcaption></figure>

J. Sutter et al. 2024). Although they do not cover exactly the same wavelengths, the WISE 3 (~12 µm) and WISE 4 (~22 µm) bands can be used for similar purposes. Additionally, because WISE performed an all-sky survey, a PAH fraction tracer constructed from these two bands would be available for all WISE-detected galaxies. To facilitate this, we examine the relation between WISE 3/WISE 4 ratios and the model parameter $q_{PAH}$.

In Figure 14 we show the relation between the ratio of WISE 3-to-WISE 4 emission and the value of $q_{PAH}$ in all of the fitted pixels of our sample (on resolved scales, above the 3σ detection). The white circles are the median values of $q_{PAH}$ in bins of $log_{10}$(WISE 3/WISE 4).

For consistency with the rest of the analysis, we perform linear fits to these medians, selecting only those with at least 250 points per bin, and show the result of that fit in the figure. This fit is described by the following equation:

$$q_{PAH} = (4.0 \pm 0.4) \times \log_{10}(\text{WISE 3/WISE 4}) + (3.8 \pm 0.2). \qquad (10)$$

However, there is clearly variation not consistent with a simple linear function visible in the binned medians.

Going further, for the sake of providing a predictive functional form only (without physical interpretation), we fit a quadratic function, shown with solid red line, to that same distribution of medians. These are marked with red diamonds. For the purpose of estimating $q_{PAH}$ with WISE bands, for values of the WISE 3/WISE 4 ratio between 0.09 and 0.9, which covers the majority of the data set, we suggest

$$q_{PAH} = (6.6 \pm 0.2) \times \log_{10}(\text{WISE 3/WISE 4})^2 + (12.7 \pm 0.3) \times \log_{10}(\text{WISE 3/WISE 4}) \times (6.2 \pm 0.1). \qquad (11)$$

## 6.5. Applications to High-redshift Studies

Observations of distant galaxies are generally limited in their spatial resolution and wavelength coverage for studying dust emission. At long wavelengths, many studies rely on ~two rest-frame far-IR measurements to fit an SED (e.g., J. Witstok et al. 2023) to constrain both dust mass and temperature for a

galaxy. At shorter wavelengths, JWST has enabled sensitive observations of mid-IR emission from galaxies out to redshifts of $z \sim 3-4$ (K. Ronayne et al. 2023; S. Alberts et al. 2024; I. Shivaei et al. 2024), but does not cover the ~20 µm continuum emission, leaving diagnostics of PAH fraction uncertain. For measuring dust mass and PAH content of high-z galaxies from their infrared emission, the key missing information is the intensity and distribution of the radiation field. Since many high-z galaxies have measurements of $M_*$ and SFR, and often ΔMS by comparison to the galaxy population at that redshift, it may be possible to use our scaling relationships to infer the radiation field properties for these galaxies. A known radiation field distribution would enable more accurate estimates of $M_d$ and $q_{PAH}$ from limited mid- or far-IR SED coverage.

A key question is whether it is possible to extrapolate the scaling relations we have determined for $z=0$ galaxies for use at higher redshifts. To test this, we would ideally need similar SED coverage for high-z galaxies, in order to apply the same type of modeling and derive radiation field properties. While most individual galaxies only have sparse SED coverage, it is possible to stack galaxies to obtain average SED constraints. M. Béthermin et al. (2015) used this approach for galaxies at $0 < z < 4$ from the COSMOS field using data from surveys of the field by Spitzer, Herschel, LABOCA, and AzTEC. They performed stacks for a sample of main-sequence galaxies at each redshift, as well as a sample of starbursts, defined to be 10 times the sSFR of the main sequence at that redshift. With the stacked SEDs, they fit B. T. Draine & A. Li (2007) models, allowing them to obtain $U_{min}$, $\gamma$, and $\bar{U}$. These results are directly comparable to our work. M. Béthermin et al. (2015) also provided average values of SFR, $M_*$, and sSFR for the stacked samples.

To test our scaling relationships and their ability to predict radiation field for high-z targets, in Figure 15 we show our integrated galaxy measurements of $\bar{U}$ as a function of SFR, $M_*$, and sSFR along with the binned values for main-sequence galaxies from $0 < z < 4$ from M. Béthermin et al. (2015). We find an excellent correspondence between their measured $\bar{U}$ and our predicted $\bar{U}$-sSFR scaling relationship, across redshift bins. Our SFR scaling does not provide a good prediction for $\bar{U}$, indicating that galaxy-averaged SFR does not trace the conditions that change the radiation field distribution. This makes sense given that a high SFR could result from either a large amount of distributed low-intensity star formation, or a small, high-intensity region of star formation. $M_*$ also fails at predicting $\bar{U}$ for the high-redshift samples, though this is also the case for the $z=0$ galaxies as well.

The correlation of dust temperature and radiation field with sSFR has been highlighted in a variety of studies in addition to M. Béthermin et al. (2015). B. Magnelli et al. (2014) used stacked Herschel observations of galaxies in deep fields to measure average dust temperature in SFR-$M_*$ space. They found that dust temperature was more strongly correlated with sSFR and ΔMS than with SFR and $M_*$. C. Schreiber et al. (2018) also fit SEDs of stacked galaxies, reproducing the trend of dust temperature increasing with redshift and showing a decreasing trend of PAH fraction (traced by their IR8 parameter) as a function of redshift and metallicity. C. Schreiber et al. (2018) also investigated trends relative to the main sequence and found correlations between the dust temperature and PAH fraction with the equivalent of ΔMS. The cause of trends of dust
<figure><img src="image_19.png" /><figcaption>**Figure 15.** Distribution of SFR, $M_*$, and sSFR as a function of the integrated fit value of $\bar{U}$, for the sample in this paper (gray points), with the best-fit linear relations shown as solid black lines (Figure 5 and Table 3). The purple symbols are from the main-sequence sample of M. Béthermin et al. (2015).</figcaption></figure>

temperature with sSFR or $\Delta$MS have been explained through combinations of increased star formation efficiency (leading to higher SFR and radiation field intensity for a given $\Sigma_d$), more compact star-forming regions, or changes in the DGR or dust-to-metals ratio leading to more pervasive, higher-intensity radiation fields for a given SFR (H. Hirashita & I.-D. Chiang 2022; L. Sommovigo et al. 2022).

The fact that local Universe sSFR scaling does a good job of predicting $\bar{U}$ up to $z \sim 4$ is surprising, considering the very different properties of galaxy disks expected at that era. Galaxy disks and star-forming clumps show high velocity dispersions and high $\Sigma_{SFR}$ along with disk scale heights that are thought to be much more vertically extended compared to $z \sim 0$ galaxies (L. J. Tacconi et al. 2020). Given two galaxy disks with similar sSFR$_R$, we might expect that the disk with the larger scale height might have lower $\bar{U}$, simply by a geometric argument. However, galaxies at higher redshift may also have overall lower metallicity and DGR, meaning that the same $\Sigma_{SFR}$ could produce a harder, more pervasive radiation field that can heat the available dust up to higher temperatures.

## 7. Conclusions

We present an analysis of dust and radiation field properties of a sample of more than 800 galaxies in the local Universe based on WISE and Herschel SED fitting. Using the B. T. Draine & A. Li (2007) dust model and the DustBFF grid-based Bayesian fitting code, we derive the properties of the distribution of radiation fields heating the dust ($U_{min}$, $\gamma$, $\bar{U}$), the dust mass surface density ($\Sigma_d$), and the fraction of dust in the form of PAHs ($q_{PAH}$) on both integrated and resolved scales. Our sample provides one of the largest, most diverse samples of galaxies where resolved maps of $\Sigma_d$ and $q_{PAH}$ are available.

We fit scaling relations between the dust model parameters with SFR, $M_*$, sSFR, and offset from the main sequence, $\Delta$MS, on integrated scales (Figures 5, 6) and resolved scales (Figures 7, 8). Results of our analysis include the following:

1. We find that all parameters describing the radiation field ($\bar{U}$, $U_{min}$, $\gamma$) are positively correlated with quantities related to star formation (SFR, sSFR, $\Delta$MS) on both integrated and resolved scales. On resolved scales, the

correlation of $U_{min}$ and $\bar{U}$ with both $\Sigma_*$ and $\Sigma_{SFR}$ suggest that the average radiation field has contributions from both old and young stellar populations.

2. The derived radiation field parameters show behavior consistent with expectations for the delta-function plus power-law model, where the high-intensity radiation fields in the vicinity of star-forming regions populate the power-law distribution, and $U_{min}$ describes the minimum interstellar radiation field. We find $\gamma$ is generally small and shows weak or no correlation with stellar mass on either resolved or integrated scales, but strong correlation with SFR related quantities.

3. We find a strong anticorrelation between the PAH fraction and sSFR and $\Delta$MS on both integrated and resolved scales. For resolved $q_{PAH}$, an increase of 1 dex in sSFR$_R$ results in $q_{PAH}$ dropping from the 4%–5% typical for solar metallicity ISM to $\sim 1\%$, comparable to what is seen in low-metallicity galaxies. Resolved $q_{PAH}$ is also strongly anticorrelated with $\Sigma_{SFR}$, emphasizing that the local star formation is a key agent affecting the PAH fraction, likely due to destruction of PAHs in H II regions. The strong anticorrelation of $q_{PAH}$ with sSFR suggests that careful sample selection is necessary to isolate trends of $q_{PAH}$ with metallicity, since many low-Z targets have been selected to highly star-forming dwarfs, and the shape of the SF main sequence already introduces a correlation of sSFR with $M_*$ (and therefore metallicity).

4. We find a constant ratio of $M_d/M_*$ with $M_*$ on integrated scales for main-sequence galaxies. Galaxies with low sSFR, falling below the main sequence, have lower $M_d/M_*$, most likely due to low gas fractions. On the main sequence, we see a constant $M_d/M_*$ value over $>2$ orders of magnitude in $M_*$, suggesting that trends of metallicity (and DGR) and gas fraction with $M_*$ cancel each other out in setting $M_d/M_*$.

5. The constant $M_d/M_*$ does not persist on resolved scales, where we observe a negative correlation of $\Sigma_d/\Sigma_*$ with $\Sigma_*$. We find that $\Sigma_d/\Sigma_*$ increases at fixed $\Sigma_*$ above the resolved main sequence. This suggests that higher $\Sigma_d/\Sigma_*$ results from higher gas (and dust) surface densities, correlated with increased $\Sigma_{SFR}$. The overall sublinear correlation of $\Sigma_d$ with $\Sigma_*$ seems to be a general property of the distribution of the ISM in galaxies, and provides an
interesting contrast to the widely observed linear CO--Σ<sub>*</sub> correlation, highlighting systematic trends in α<sub>CO</sub>.

6. We find that Σ<sub>SFR</sub> ∝ Σ<sub>d</sub><sup>1.4±0.08</sup>. The slope of ~1.4 suggest our Σ<sub>d</sub> measurements extend into the H I-dominated parts of galaxies. We measure a typical “dust depletion time,” analogous to the more standard gas or H<sub>2</sub> depletion time, of ≈8.7 Myr. This value of SFR per unit dust mass may be a useful comparison for high-z systems where gas masses, particularly atomic gas, may be more challenging to measure than dust masses.

7. Using a combination of Σ<sub>SFR</sub> and Σ<sub>*</sub>, our 2D scaling relationships (Section 5.6) are able to predict <U> to an RMSE of 0.2 dex, q<sub>PAH</sub> to 1.2%, and Σ<sub>d</sub> to 0.26 dex.

8. We assess the predictive power of our integrated-galaxy radiation field scaling relations for high-redshift galaxies, using stacked SEDs of main-sequence galaxies from M. Béthermin et al. (2015). We find that our <U>-sSFR scaling relationship matches the <U> values for the stacked main-sequence galaxies from 0 < z < 4.

Associated with this paper, we deliver all reduced Herschel data and dust fitting results to the Infrared Science Archive (IRSA), as described in Appendix A (doi:10.26131/IRSA581). We note that our delivery includes Herschel data for a substantial number of galaxies that were not included in our analysis, due to lacking both PACS and SPIRE coverage (or WISE data in the z0MGS delivery), and it includes the SPIRE 350 and 500 µm bands that we did not use in our SED fitting. We also provide the resolution-matched, background-subtracted maps for all galaxies for which we perform SED fitting. These data products can be straightforwardly ingested into other SED-fitting procedures that make different assumptions about the dust and radiation field properties.

Until the next far-infrared space telescope is launched, the set of nearby galaxies observed by Herschel, that we compile here, will remain a uniquely important resource for studying interstellar dust in galaxies.

## Acknowledgments

We thank the referee for providing a careful reading of the paper and suggestions that improved its clarity. The work of J. C., K.S., I.C., and A.K.L. is supported by NASA ADAP grants NNX16AF48G and NNX17AF39G and National Science Foundation grant No. 1615728. J.C. acknowledges support from ERC starting grant No. 851622 DustOrigin, and funding from the Belgian Science Policy Office (BELSPO) through the PRODEX project “JWST/MIRI Science exploitation” (C4000142239). I.C. acknowledges the support from the National Science and Technology Council grant 111-2112-M-001-038-MY3, and the Academia Sinica investigator award AS-IA-109-M02 (PI: Hiroyuki Hirashita). E.W.K. acknowledges support from the Smithsonian Institution as a Sub-millimeter Array (SMA) Fellow and the Natural Sciences and Engineering Research Council of Canada.

We are grateful to the staff of IPAC for their support in hosting the large data set presented in this paper. We acknowledge the usage of the HyperLeda database (http://leda.univ-lyon1.fr). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France (M. Wenger et al. 2000), TOPCAT, an interactive graphical viewer and editor for tabular data (M. B. Taylor 2005), pandas (W. McKinney 2010, 2011), APLpy, an open-source plotting

package for Python hosted at http://aplpy.github.io, SciPy (P. Virtanen et al. 2020), ds9, a tool for data visualization supported by the Chandra X-ray Science Center (CXC) and the High Energy Astrophysics Science Archive Center (HEA-SARC) with support from the JWST Mission office at the Space Telescope Science Institute for 3D visualization, matplotlib, a Python library for publication quality graphics (J. D. Hunter 2007), Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018, 2022), and NumPy (C. R. Harris et al. 2020). This publication makes use of data products from the Wide-field Infrared Survey Explorer (E. L. Wright et al. 2010), which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including Univ. Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sussex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KUL, CSL, IMEC (Belgium); CEA, OAMP (France); MPIA (Germany); IFSI, OAP/AOT, OAA/CAISMI, LENS, SISSA (Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/ CNES (France), DLR (Germany), ASI (Italy), and CICT/MCT (Spain). This research has made use of HIPE, a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS, and SPIRE consortia.

Facilities: WISE, Herschel

## Appendix A
Data Delivery

### A.1. Catalog

We provide a machine-readable table containing the parameters describing each galaxy, the integrated photometry, and our integrated SED-fitting results for each target. These are described in Table 1 and in detail, this includes: galaxy name; PGC name; R.A. and decl.; *r*<sub>25</sub> value; effective radius *r*<sub>e</sub>; inclination; position angle, distance, integrated flux in bands from WISE 1 to SPIRE 250 and associated errors; the type of fit available (resolved or integrated, where resolved means integrated is available as well); and the SED-fitting parameter results including *U*<sub>min</sub>, <U>, γ, *q*<sub>PAH</sub>, and *M*<sub>d</sub>, and associated errors.

### A.2. Resolved Maps

We provide four sets of FITS files related to the resolved analysis from this work:

1. [GalaxyName]_[Band].fits: the Herschel images after Scanamorphos reduction and the background subtraction at native resolution, in MJy sr⁻¹. The header of these maps contains Scanamorphos keywords, and includes the observation IDs from the Herschel archive used to create the
final map; In the header of these maps, we add a few key elements that were used in the data processing:

(a) R25COEFF: the A coefficient used to scale $r_{25}$ to create the galaxy mask (Section 3.3);

(b) COEFF[1,2,3]: the coefficients of the plane to remove a 2D background;

(c) ADDGAL[x]: the name(s) of any galaxy that was found in the cutout and consequently masked following identical procedure as described in Section 3.3;

2. [GalaxyName]_[Band]_conv250.fits: the background-subtracted, convolved to SPIRE 250 18" resolution, regridded Herschel maps, in MJy sr$^{-1}$. The headers of these maps contain the same extra keywords as the ones mentioned just above;

3. [GalaxyName]_DustParameters Conv250. fits: a multi-extensions file containing the realizations maps (K. D. Gordon et al. 2014) of the B. T. Draine & A. Li (2007) dust parameters, $U_{min}$, $\gamma$, $q_{PAH}$, $\Sigma_d$, and $\Omega_*$, as well as $\bar{U}$ (not fitted, calculated), for pixels passing a 1$\sigma$ S/N cut in the data (Section 4.2), and their associated

16th–84th percentiles. The maps and errors are calculated from 100 realization maps for each galaxy;

4. [GalaxyName]_extra.fits: a multi-extension file containing:

(a) the galaxy mask, with the used R25COEFF coefficient;

(b) masks of the pixels passing the $1\sigma$, $2\sigma$, and $3\sigma$ S/N cuts;

(c) the final “master-mask” at the SPIRE 250 resolution, masking the galaxy, stars, and bright sources.

## Appendix B: Dust Emission Parameters Corner Plot

Figure 16 shows the 2D and 1D marginalized histograms of the fitted dust and radiation field parameters and their relationships with each other. The clear relation between $U_{min}$ and $\bar{U}$ is expected due to the definition of $\bar{U}$, see Equation (7). The scatter in the $\bar{U}-U_{min}$ plot is correlated with $\gamma$. There is a degree of expected degeneracy between the model parameters

<figure><img src="image_20.png" /><figcaption>Figure 16. Distributions of the fitted dust and radiation field model parameters for all resolved pixel in this analysis passing the cuts (Section 4.4). The top plot of each column shows a marginalized 1D histogram of that parameter. The other panels show 2D histograms of all pixels in the resolved analysis. Note that $\bar{U}$ is not a fitted parameter, but is derived from the combination of $U_{min}$ and $\gamma$.</figcaption></figure>
(extensively investigated in F. Galliano et al. 2021), which is visible in the Figure: the far-IR SED responds strongly to a change in $U_{min}$, shifting the IR peak to shorter wavelengths, and in $\Sigma_d$, scaling the SED vertically. This is seen in the anticorrelation between these two parameters, and is similar to the $\beta-T_d$ relation investigated in the literature (B. C. Kelly et al. 2012; F. Galliano et al. 2021). Both of these parameters, and the stellar light scaling factor, $\Omega_*$, also affect the mid-IR SED, and influence the value of $q_{PAH}$.

The diagonal cut visible in the $U_{min}-\Sigma_d$ and $\bar{U}-\Sigma_d$ panels is due to selection effects. The typical sensitivity of the Herschel maps can be translated into a combination of $\bar{U}$ and $\Sigma_d$. For low $\bar{U}$, the only detectable dust emission will correspond to high $\Sigma_d$ values. At high $\bar{U}$, even small $\Sigma_d$ can potentially be detected. This results in the diagonal cut clearly seen in the $\bar{U}-\Sigma_d$ and $U_{min}-\Sigma_d$ panels.

The $\Omega_*$ parameter scales a 5000 K blackbody to account for starlight, primarily affecting the shortest mid-IR wavelengths. For the purposes of our fitting, we chose to use a coarse grid in $\Omega_*$ since it was not critical to our analysis. The coarseness of the $\Omega_*$ grid allowed us to use finer sampling of the other SED parameters. The visible stripes in the bottom row of panels are due to this coarse sampling.

## Appendix C
### Effect of the Wavelength Coverage

As mentioned in the text, we only use data up to SPIRE 250, as a compromise between spectral coverage and resolution. In the range of temperature (radiation field) covered in this sample, the infrared peak will fall short of 250 µm, and using SPIRE 250 only is sufficient to measure the far-IR peak. Several works have shown the effect of wavelength coverage and resolution on the recovery of dust parameters (G. Aniano et al. 2012). The effect of spatial resolution was also investigated by F. Galliano et al. (2011).

As a test of the accuracy of fitting up to SPIRE 250, Figure 17 shows the distributions of the dust parameters for two fits at the same resolution of SPIRE 500, ~36", using different spectral coverage: up to 500 µm and up to 250 µm. Data and results are those from J. Chastenet et al. (2021), covering the mid- to far-IR emission of M101 (NGC 5457). The top row shows the pixel-by-pixel values for the two different fits. Only the radiation field values show a clear deviation from the 1:1 line, and it becomes significant (slope of ~1.25) only at high $U_{min}$ values. The bottom rows show that the distributions of values in both cases are virtually identical.

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Distributions of the dust parameter values for M101 using the data from J. Chastenet et al. (2021). These values are extracted after fitting mid- to far-IR data at the SPIRE 500 resolution (~36") up to 500 µm, or up to 250 µm only. Top row: values from a fit using all data up to 500 µm (x-axis) vs. values from a fit using data up to 250 µm only (y-axis). The black dashed line shows the 1:1 relation. Only the radiation field is noticeably affected by wavelength coverage, and the difference becomes large only at high $U_{min}$ values. Bottom row: histograms for the dust parameters.</figcaption></figure>
ORCID iDs

Jérémy Chastenet https://orcid.org/0000-0002-5235-5589

Karin Sandstrom https://orcid.org/0000-0002-4378-8534

Adam K. Leroy https://orcid.org/0000-0002-2545-1700

Caroline Bot https://orcid.org/0000-0001-6118-2985

I-Da Chiang (江宜達) https://orcid.org/0000-0003-2551-7148

Ryan Chown https://orcid.org/0000-0001-8241-7704

Karl D. Gordon https://orcid.org/0000-0001-5340-6774

Eric W. Koch https://orcid.org/0000-0001-9605-780X

Jessica Sutter https://orcid.org/0000-0002-9183-8102

Thomas G. Williams https://orcid.org/0000-0002-0012-2142

References

Abdurro'uf, Lin, Y.-T., Hirashita, H., et al. 2022, ApJ, 935, 98

Alberts, S., Lyu, J., Shivaei, I., et al. 2024, ApJ, 976, 224

Allamandola, L. J., Tielen, A. G. G. M., & Barker, J. R. 1989, ApJS, 71, 733

Aniano, G., Draine, B. T., Calzetti, D., et al. 2012, ApJ, 756, 138

Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218

Aniano, G., Draine, B. T., Hunt, L. K., et al. 2020, ApJ, 889, 150

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Baker, W. M., Maiolino, R., Bluck, A. F. L., et al. 2022, MNRAS, 510, 3622

Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Bianchi, S., De Vis, P., Viaene, S., et al. 2018, A&A, 620, A112

Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

Bigiel, F., Leroy, A., Walter, F., et al. 2010, AJ, 140, 1194

Blanc, G. A., Heiderman, A., Gebhardt, K., Evans, N. J., II, & Adams, J. 2009, ApJ, 704, 842

Blitz, L., & Rosolowsky, E. 2006, ApJ, 650, 933

Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207

Boquien, M., Buat, V., Boselli, A., et al. 2012, A&A, 539, A145

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Boquien, M., & Salim, S. 2021, A&A, 653, A149

Boselli, A., Eales, S., Cortese, L., et al. 2010, PASP, 122, 261

Calura, F., Pozzi, F., Cresci, G., et al. 2017, MNRAS, 465, 54

Casasola, V., Bianchi, S., De Vis, P., et al. 2020, A&A, 633, A100

Casasola, V., Bianchi, S., Magrini, L., et al. 2022, A&A, 668, A130

Casasola, V., Cassarà, L. P., Bianchi, S., et al. 2017, A&A, 605, A18

Casasola, V., Hunt, L., Combes, F., & García-Burillo, S. 2015, A&A, 577, A135

Cesarsky, D., Lequeux, J., Abergel, A., et al. 1996, A&A, 315, L309

Chastenet, J., Sandstrom, K., Chiang, I. D., et al. 2019, ApJ, 876, 62

Chastenet, J., Sandstrom, K., Chiang, I.-D., et al. 2021, ApJ, 912, 103

Chastenet, J., Sutter, J., Sandstrom, K., et al. 2023, ApJL, 944, L11

Chiang, I.-D., Hirashita, H., Chastenet, J., et al. 2023, MNRAS, 520, 5506

Chiang, I.-D., Sandstrom, K. M., Chastenet, J., et al. 2024, ApJ, 964, 18

Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128

Clark, C. J. R., De Vis, P., Baes, M., et al. 2019, MNRAS, 489, 5256

Clark, C. J. R., Dunne, L., Gomez, H. L., et al. 2015, MNRAS, 452, 397

Clark, C. J. R., Roman-Duval, J. C., Gordon, K. D., Bot, C., & Smith, M. W. L. 2021, ApJ, 921, 35

Clark, C. J. R., Roman-Duval, J. C., Gordon, K. D., et al. 2023, ApJ, 946, 42

Clark, C. J. R., Verstocken, S., Bianchi, S., et al. 2018, A&A, 609, A37

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

Cortese, L., Ciesla, L., Boselli, A., et al. 2012, A&A, 540, A52

Croxall, K. V., Smith, J. D., Wolfire, M. G., et al. 2012, ApJ, 747, 81

da Cunha, E., Eminian, C., Charlot, S., & Blaizot, J. 2010, MNRAS, 403, 1894

Dale, D. A., Boquien, M., Turner, J. A., et al. 2023, AJ, 165, 260

Dale, D. A., Cohen, S. A., Johnson, L. C., et al. 2009, ApJ, 703, 517

Dale, D. A., & Helou, G. 2002, ApJ, 576, 159

Dale, D. A., Helou, G., Contursi, A., Silbermann, N. A., & Kolhatkar, S. 2001, ApJ, 549, 215

Davies, J. I., Baes, M., Bendo, G. J., et al. 2010, A&A, 518, L48

Davies, J. I., Baes, M., Bianchi, S., et al. 2017, PASP, 129, 044102

Davies, J. I., Nersesian, A., Baes, M., et al. 2019, A&A, 626, A63

De Looze, I., Fritz, J., Baes, M., et al. 2014, A&A, 571, A69

De Looze, I., Lamperti, I., Saintonge, A., et al. 2020, MNRAS, 496, 3668

De Vis, P., Dunne, L., Maddox, S., et al. 2017a, MNRAS, 464, 4680

De Vis, P., Gomez, H. L., Schofield, S. P., et al. 2017b, MNRAS, 471, 1743

De Vis, P., Jones, A., Viaene, S., et al. 2019, A&A, 623, A5

Decleir, M., De Looze, I., Boquien, M., et al. 2019, MNRAS, 486, 743

Dowell, C. D., Pohlen, M., Pearson, C., et al. 2010, Proc. SPIE, 7731, 773136

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T., & Anderson, N. 1985, ApJ, 292, 494

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Draine, B. T., Li, A., Hensley, B. S., et al. 2021, ApJ, 917, 3

Egorov, O. V., Kreckel, K., Sandstrom, K. M., et al. 2023, ApJL, 944, L16

Eibensteiner, C., Sun, J., Bigiel, F., et al. 2024, A&A, 691, A163

Ellison, S. L., Sánchez, S. F., Ibarra-Medel, H., et al. 2018, MNRAS, 474, 2039

Engelbracht, C. W., Gordon, K. D., Rieke, G. H., et al. 2005, ApJL, 628, L29

Engelbracht, C. W., Rieke, G. H., Gordon, K. D., et al. 2008, ApJ, 678, 804

Galliano, F. 2022, Habilitation thesis, Université Paris-Saclay

Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673

Galliano, F., Hony, S., Bernard, J. P., et al. 2011, A&A, 536, A88

Galliano, F., Nersesian, A., Bianchi, S., et al. 2021, A&A, 649, A18

Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69

Gordon, K. D., Engelbracht, C. W., Rieke, G. H., et al. 2008, ApJ, 682, 336

Gordon, K. D., Roman-Duval, J., Bot, C., et al. 2014, ApJ, 797, 85

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Helou, G. 1986, ApJL, 311, L33

Hensley, B. S., & Draine, B. T. 2021, ApJ, 906, 73

Hensley, B. S., & Draine, B. T. 2023, ApJ, 948, 55

HeVICS Team 2020, Herschel Virgo Cluster Survey, IPAC, doi:10.26131/IRSA70

Hirashita, H., & Chiang, I.-D. 2022, MNRAS, 516, 1612

Hsieh, B. C., Lin, L., Lin, J. H., et al. 2017, ApJL, 851, L24

Hunt, L. K., De Looze, I., Boquien, M., et al. 2019, A&A, 621, A51

Hunt, L. K., Draine, B. T., Bianchi, S., et al. 2015, A&A, 576, A33

Hunter, J. D. 2007, CSE, 9, 90

Jenkins, E. B. 2009, ApJ, 700, 1299

Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46

Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, ApJ, 752, 55

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Kennicutt, R. C., Calzetti, D., Aniano, G., et al. 2011, PASP, 123, 1347

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Law, K.-H., Gordon, K. D., & Misselt, K. A. 2018, ApJS, 236, 32

Law, K.-H., Gordon, K. D., & Misselt, K. A. 2021, ApJ, 920, 96

Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170

Leroy, A. K., Sandstrom, K. M., Lang, D., et al. 2019, ApJS, 244, 24

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782

Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013, AJ, 146, 19

Li, A. 2020, NatAs, 4, 339

Lin, L., Pan, H.-A., Ellison, S. L., et al. 2019, ApJL, 884, L33

Madden, S. C., Galliano, F., Jones, A. P., & Sauvage, M. 2006, A&A, 446, 877

Madden, S. C., Rémy-Ruyer, A., Galametz, M., et al. 2013, PASP, 125, 600

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Maiolino, R., & Mannucci, F. 2019, A&ARv, 27, 3

Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13

Marble, A. R., Engelbracht, C. W., van Zee, L., et al. 2010, ApJ, 715, 506

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212

McKinney, W. 2010, Proc. 9th Python in Science Conf., ed. S. van der Walt & J. Millman, (SciPy) 51

McKinney, W. 2011, Python for High Performance and Scientific Computing 14 https://www.researchgate.net/publication/265194455_pandas_a_Foundational_Python_Library_for_Data_Analysis_and_Techniques

Muñoz-Mateos, J. C., Gil de Paz, A., Boissier, S., et al. 2009, ApJ, 701, 1965

Nersesian, A., Verstocken, S., Viaene, S., et al. 2020, A&A, 637, A25

Nersesian, A., Xilouris, E. M., Bianchi, S., et al. 2019, A&A, 624, A80

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Ott, S. 2010, in ASP Conf. Ser. 434, Astronomical Data Analysis Software and Systems XIX, ed. Y. Mizumoto, K.-I. Morita, & M. Ohishi (San Francisco, CA: ASP) 139
Paladini, R., Umana, G., Veneziani, M., et al. 2012, ApJ, 760, 149

Paradis, D., Mény, C., Demyk, K., Noriega-Crespo, A., & Ristorcelli, I. 2023, A&A, 674, A141

Paradis, D., Paladini, R., Noriega-Crespo, A., et al. 2011, ApJ, 735, 6

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Pessa, I., Schinnerer, E., Leroy, A. K., et al. 2022, A&A, 663, A61

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A132

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Povich, M. S., Stone, J. M., Churchwell, E., et al. 2007, ApJ, 660, 346

Relaño, M., Kennicutt, R., Lisenfeld, U., et al. 2016, A&A, 595, A43

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2014, A&A, 563, A31

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121

Roman-Duval, J., Jenkins, E. B., Tchernyshyov, K., et al. 2022, ApJ, 935, 105

Ronayne, K., Papovich, C., Yang, G., et al. 2023, ApJ, 970, 61

Roussel, H. 2013, PASP, 125, 1126

Rushton, M. T., Popescu, C. C., Inman, C., Natale, G., & Pricopi, D. 2022, MNRAS, 514, 113

Saikia, G., Patra, N. N., Roy, N., & Jog, C. J. 2020, MNRAS, 492, 2517

Saintonge, A., & Catinella, B. 2022, ARA&A, 60, 319

Saintonge, A., Wilson, C. D., Xiao, T., et al. 2018, MNRAS, 481, 3497

Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529

Sánchez, S. F. 2020, ARA&A, 58, 99

Sánchez, S. F., Rosales-Ortega, F. F., Iglesias-Páramo, J., et al. 2014, A&A, 563, A49

Sandstrom, K. M., Bolatto, A. D., Draine, B. T., Bot, C., & Stanimirović, S. 2010, ApJ, 715, 701

Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30

Schruba, A., Leroy, A. K., Walter, F., et al. 2011, AJ, 142, 37

Schruba, A., Leroy, A. K., Walter, F., et al. 2012, AJ, 143, 138

Sellgren, K. 1984, ApJ, 277, 623

Shetty, R., Kelly, B. C., Rahman, N., et al. 2014, MNRAS, 437, L61

Shivaei, I., Alberts, S., Florian, M., et al. 2024, A&A, 690, A89

Smith, J. D. T., Draine, B. T., Dale, D. A., et al. 2007, ApJ, 656, 770

Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2020, ApJ, 892, 148

Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2023, ApJL, 945, L19

Sutter, J., Sandstrom, K., Chastenet, J., et al. 2024, ApJ, 971, 178

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Taylor, M. B. 2005, in ASP Conf. Ser. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert (San Francisco, CA: ASP), 29

Teng, Y.-H., Chiang, I.-D., Sandstrom, K. M., et al. 2024, ApJ, 961, 42

Tielens, A. G. G. M. 2008, ARA&A, 46, 289

Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898

Utomo, D., Chiang, I.-D., Leroy, A. K., Sandstrom, K. M., & Chastenet, J. 2019, ApJ, 874, 141

Viaene, S., Baes, M., Tamm, A., et al. 2017, A&A, 599, A64

Viaene, S., Fritz, J., Baes, M., et al. 2014, A&A, 567, A71

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vutisalchavakul, N., Evans, N. J., II, & Battersby, C. 2014, ApJ, 797, 77

Wang, J., Koribalski, B. S., Serra, P., et al. 2016, MNRAS, 460, 2143

Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9

Whitcomb, C. M., Smith, J. D. T., Sandstrom, K., et al. 2024, ApJ, 974, 20

Wieprecht, E., Schreiber, J., de Jong, J., et al. 2009, in ASP Conf. Ser. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler (San Francisco, CA: ASP), 531

Witstok, J., Jones, G. C., Maiolino, R., Smit, R., & Schneider, R. 2023, MNRAS, 523, 3119

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Ysard, N., Jones, A. P., Guillet, V., et al. 2024, A&A, 684, A34

z0MGS Team 2019, z=0 Multiwavelength Galaxy Synthesis, IPAC, doi:10.26131/IRSA6