The Fraction of Dust Mass in the Form of PAHs on 10–50 pc Scales in Nearby Galaxies

This work uses data from three of the MIRI bands: F770W, F1130W, and F2100W (Rieke et al., 2015) as well as the NIRCam F200W. The MIRI F770W and F1130W filters capture the emission from two of the strongest PAH features, the 7.7 µm and 11.3 µm features. The F2100W filter is dominated by dust continuum emission. As an example, the F770W and F2100W maps of NGC 4303 are shown in Figure 1. The NIRCam F200W filter primarily covers stellar emission and is therefore used to remove starlight from the F770W data (see Section 3.1).

Each target was covered by small 1- to 4-pointing mosaics, chosen to maximize the overlap with the available optical and millimeter spectroscopy and cover the majority ($\geq 70\%$) of star formation activity traced by MIR emission in each galaxy (defined by a contour of 0.5 MJy sr^-1 at WISE 12 µm). Details of the data processing are described in Williams et al. (2024). Briefly, JWST data are reduced using the most up-to-date calibration reference data system (CRDS) context (jwst_1201), and JWST calibration pipeline version v1.13.4, and the PHANGS-JWST data version 1.1.0. In addition to the calibration pipeline, the background level in each filter is matched to existing wide-field Spitzer or WISE data that extend off of the galaxy, as discussed in Leroy et al. (2023b). The JWST data are finally all convolved to a common Gaussian resolution of 0$\farcs$90, a resolution slightly larger than the F2100W beam, to improve signal-to-noise at F2100W and match the resolution of existing multiwavelength data (see discussion in Williams et al., 2024). This corresponds to only a moderate decrease in resolution (native resolution of 0$\farcs$67 at F2100W, the largest of this filter set) but greatly improves the signal-to-noise, as the noise in these filters is pixel based rather than resolution element based. The convolution was implemented using the method described in Aniano et al. (2011), using the WebbPSF²²2https://stsci.app.box.com/v/jwst-simulated-psf-library models to generate kernels³³3https://github.com/francbelf/jwst_kernels. The pixel scale is left at the original pixel scale of 0$\farcs$11 which over samples the PSF, but this does not affect the analysis presented in this paper.

A 3$\sigma$ signal-to-noise cut is performed on each MIRI filter used in this work. We determine the typical noise by finding empty sky regions in a subset of the maps that were large enough to include areas off the galaxy. This was possible for: NGC 1087, NGC 1385, NGC 1433, NGC 1512, NGC 1566, and NGC 7496. The 1$\sigma$ level was determined as the average of the standard deviations of the sky regions found in these six maps. All background regions showed similar standard deviations and median values. With this method, our 3$\sigma$ limits are 0.09 MJy/sr in F770W, 0.13 MJy/sr in F1130W, and 0.29 MJy/sr in F2100W. These values are in good agreement with the characterization of noise levels presented in Williams et al. (2024) and match well with the expected values from convolving the pipeline-generated error maps. Additionally, several of the nuclei in our sample saturated in F2100W, producing diffraction spikes across the maps. By-eye masking is done to exclude regions where these spikes contaminate the maps. In convolving to lower resolution, we also convolve the spike masks and use them to conservatively exclude any contaminated pixels by removing all pixels where the convolved mask has a value greater than 0.9.

In order to assess the combination of JWST bands as tracers of dust properties, we compare the JWST data to the results of MIR to FIR spectral energy distribution fitting to the Draine & Li (2007) dust models produced as part of the “$z=0$ Multiwavelength Galaxy Synthesis” program (z0MGs, Leroy et al., 2019; Chastenet et al., 2021, Chastenet et al, in prep). The results of the Draine & Li (2007) model fitting include maps of $q_{\rm{PAH}}$, which we use to calibrate our JWST photometric measurement of the PAH fraction, $R_{\rm PAH}$. These fit results were produced using a combination of WISE and Herschel Space Observatory data and the Draine & Li (2007) dust models (using the DustBFF code, Gordon et al., 2014).

In addition to the infrared data, we use ancillary data products to determine the gas and stellar properties of the galaxies in our sample. These include ${}^{12}{\rm CO}~(2-1)$ moment 0 maps from PHANGS–ALMA (Leroy et al., 2021b, a), H$\alpha$ maps from PHANGS–MUSE (Emsellem et al., 2022), and HI maps from several VLA projects including THINGS (Walter et al., 2008) or MeerKAT (Sun et al., 2022), depending on data availability. MeerKAT HI maps are used in instances where both VLA and MeerKAT data are available.

The sources of HI maps are listed in Table 1. These maps have resolutions ranging from 11″ to 15″, and are left at their native resolutions for our analysis. While the typical HI resolution is much larger than the other maps used in this work, the smoothness and flatness of the atomic gas profile observed in the Milky Way and other galaxies suggests that convolution to matched resolutions will not greatly impact our comparisons (Schruba et al., 2011; Bigiel & Blitz, 2012; Leroy et al., 2013; Wong et al., 2013). The HI intensities are converted to atomic gas densities including helium and assuming optically thin HI gas, using the equation $\Sigma_{\rm{HI}}[\rm{M}_{\odot}\rm{pc}^{-2}]=0.020~I_{\rm{HI}}$[K km s^-1].

The ${}^{12}{\rm CO}~(2-1)$ moment 0 maps were observed with ALMA using the 12m, 7m, and total power array, producing maps with $\sim$1″ resolution (Leroy et al., 2021a). The ${}^{12}{\rm CO}~(2-1)$ intensities were converted to a molecular gas surface density using a constant $\alpha_{\rm CO,1-0}$ of 4.35 M_⊙ pc^-2 (K km s^-1)^-1 and a ${}^{12}{\rm CO}~(2-1)$/ ¹²CO($1-0$) line ratio $R_{21}=0.65$ (den Brok et al., 2022; Leroy et al., 2021a). For this work we require only a coarse estimate of the molecular gas surface density and so neglect conversion factor variations due to metallicity effects, excitation variations, and opacity variations. We use the broad moment 0 masks, which maximize the ${}^{12}{\rm CO}~(2-1)$ emission included in the moment 0 maps and are further described in Leroy et al. (2021a).

H$\alpha$ maps were obtained from the PHANGS–MUSE data (Emsellem et al., 2022) for all 19 galaxies. For this work, we use the convolved-optimized (“copt”) resolution maps. The copt resolution provides a uniform PSF for all MUSE data for a single galaxy using the broadest PSF for all observations. Resolutions for each MUSE map are listed in Table 1. The H$\alpha$ line was fit using the pPXF tool to fit the stellar continuum (including H$\alpha$ absorption) and a Gaussian profile for each spaxel. Further discussion of the MUSE data products can be found in Emsellem et al. (2022). In general, the MUSE “copt” resolutions are very similar to or slightly larger than the 0$\farcs$9 resolution of the convolved F2100W data. We do not do any additional steps of resolution matching, since we primarily rely on the nebular catalog for our analysis, as described below.

In addition to the maps described above, we rely on several of the PHANGS higher order data products to distinguish between environments and compare to physical properties in individual regions. These data products are described below.

H II regions are identified based on the catalogs compiled by Santoro et al. (2022) and further described in Groves et al. (2023). These nebular region catalogs were produced using the H$\alpha$ maps from the PHANGS–MUSE survey (described above) and the HIIPhot algorithm for segmentation described in Thilker et al. (2000). Below we distinguish between $R_{\rm PAH}$ inside ($R_{\rm PAH,Neb}$) and outside ($R_{\rm PAH,Diffuse}$) of these nebular regions. When we do so, we include all nebular regions in the catalogs, which are mostly H II regions but also include some other nebular regions like supernova remnants.

Metallicity maps for each galaxy were produced in Williams et al. (2022). These maps were created with the MUSE data described above (Emsellem et al., 2022) and the Pilyugin & Grebel (2016) S calibration, which relies on the relative strengths of the H$\beta$, [OIII]$\lambda 5007$, [NII]$\lambda 6584$, and [SII]$\lambda 6717+6731$ emission lines. All lines were corrected for dust extinction using the Balmer decrement before the metallicity was calculated. The Williams et al. (2022) metallicity maps interpolate between H II region metallicities using Gaussian Process Regression, yielding a fully sampled metallicity map for each galaxy.

We also use the environmental maps from Querejeta et al. (2021) to isolate morphological environments within each galaxy. These maps were created with the Spitzer 3.6µm images and primarily distinguish between centers, bars, spiral arms, inter-arm disks, and disks without strong spiral arms.

Finally, we use the extinction-corrected H$\alpha$ star formation rate (SFR) and stellar mass (M_⋆) maps produced by Belfiore et al. (2023). These maps were produced using the integral field spectroscopy from the PHANGS–MUSE survey (Emsellem et al., 2022). The H$\alpha$ measurements are corrected for attenuation using an $E(B-V)$ determined based on the measured Balmer decrement across the maps. Stellar masses were determined using Voronoi binning of the MUSE data cubes and full-spectral fitting using the pPXF tool (Cappellari & Emsellem, 2004), as described in Emsellem et al. (2022), Section 5.2.4.

To test trends in the PAH distribution, we examine the MIRI data from each galaxy on three scales: pixel-by-pixel, in circular regions with radii of 1 kpc spaced using a hexagonal grid to evenly sample the full MIRI maps, and the average of each band across the full maps. The pixel-by-pixel trends leverage the resolution of JWST and are the finest statistical sample available with this dataset. Because we are also interested in how large-scale environment affects $R_{\rm PAH}$ (e.g. $\Sigma_{\star}$, $\Sigma_{\rm{SFR}}$, sSFR) we also use larger (1 kpc) regions. Within each 1 kpc region, we record the values of $R_{\rm PAH}$ in the nebular regions ($R_{\rm PAH,Neb}$) and also diffuse gas ($R_{\rm PAH,Diffuse}$) and also note that this cannot be done pixel-by-pixel, as each pixel can only be in either a nebular region or a diffuse region. Finally, in using the full maps, we determine how integrated galaxy properties might enhance or decrease the PAH fraction. These integrated diagnostics also provide the most direct comparison to higher redshift measurements of PAHs, which will primarily yield galaxy–integrated PAH data.

In rare cases in the PHANGS sample, where there is low dust column relative to the stellar surface density (e.g. in some galactic bulges) and/or a very low PAH fraction, the F770W band can have a significant contribution from starlight. F1130W is at a long enough wavelength where the contribution from starlight will always be minimal. To demonstrate this, we show two SED models produced using the Code for Investigating GALaxy Emission (CIGALE, Boquien et al., 2019) in Figure 4. In one model, the PAH fraction is set to a relatively high value of $q_{\rm PAH}$ = 6.63% (solid lines), while in the other model the PAH fraction is set to a low value of $q_{\rm PAH}$ = 0.47% (dashed lines). The red and black lines represent the dust and total emission, respectively. Both models assume the same stellar population, with the unattenuated starlight from this population shown as the solid purple line. In the case of the low-PAH fraction model (or in the case of a model with the same $q_{\rm PAH}$ but a far lower dust column density), it is clear that the F770W MIRI filter will have a substantial contribution from starlight.

In order to remove contributions from starlight from the F770W photometry, we scale the NIRCam F200W maps to predict the amount of starlight in the F770W band, following similar procedure to what was done by Helou et al. (2004) using Spitzer 3.6 µm (see also Draine et al., 2007; Dale et al., 2009). We use a ratio of F770W_starlight/F200W = 0.13$\pm$0.04, determined from a range of reasonable stellar emission models with ages spanning 0.5 Gyr–13 Gyr, a later burst of star-formation, and the Chabrier (2003) initial mass function. In general, the results are fairly insensitive to the details of the stellar populations (see also Appendix B of Ciesla et al., 2014, for a detailed investigation). These starlight models are also assumed to experience attenuation following Calzetti et al. (2000). Dust emission is set using the Draine et al. (2014) models, with the available range of $q_{\rm PAH}$. CIGALE does not couple the attenuation curve and $q_{\rm PAH}$ assumption. It assumes an energy balance, so all light attenuated in the optical and UV is re–emitted in the infrared. CIGALE does not adjust the 2175$\AA$ feature when $q_{\rm PAH}$ changes, but for energy balance to be maintained, the integral of the infrared SED will impact the level of attenuation applied to the starlight. The stellar emission models are produced using CIGALE (Boquien et al., 2019), with all additional properties held constant using the default CIGALE settings. Synthetic photometry is performed on each stellar emission model to determine the relative fluxes in the F770W and the attenuated stellar emission model to determine the flux in the F200W bands. We use the F200W band to correct for starlight instead of the F300M due to the higher signal-to-noise of the F200W measurements and the narrower range of F770W/F200W ratios produced by different stellar emission models, as shown in Figure 5.

One potential issue with using the F200W band to correct for starlight is the possible presence of Paschen $\alpha$ emission in this filter. To confirm the small influence of the Paschen $\alpha$ line, we once again use the grid of CIGALE models produced to estimate the contribution of Paschen $\alpha$ emission to the F200W band. As CIGALE also models the nebular emission, we are able to estimate the fraction of the F200W flux from Paschen $\alpha$ emission across all of our models. Within the range of models we test, we find the maximum contribution Paschen $\alpha$ is 2%, with a median contribution of 0.04%, confirming that this emission line will not significantly affect our starlight correction using the F200W filter. It should be noted that the models tested here do not include a wide range of nebular gas emission properties, which is likely limiting the range of contributions from Paschen $\alpha$ to the F200W filter in these models. We still advocate for the use of the F200W filter as an indicator of starlight in F770W for the galaxies included in our sample due to the higher signal to noise in F200W compared to F300M.

As an additional test, star subtractions in F770W using both F200W and F300M were compared. We find that the average differences between these measurements are $\leq 0.01$ MJy/Sr, below the SNR threshold for the F770W data. This suggests that the choice of starlight filter will not impact our overall results. There are some instances where this may no longer be the case, specifically when the interstellar radiation field is dominated by an older stellar population. In the top panel of Figure 5, several repeating peaks can be seen in the F770W_⋆/F300M plot. These peaks represent models with increasingly older stellar populations, which shifts the amount of modeled starlight in the F300M band relative to the F770W band. Because the median does not reflect these older stellar populations, we find that the nuclei of our galaxies, which will be dominated by an older stellar population, starlight subtractions using F300M remove much less flux than the F200W starlight subtraction. This is a further reason for our choice of the F200W band for the starlight subtraction, which is less dependant on the age of the stellar population.

As shown above, the most obvious factor driving changes in $R_{\rm PAH}$ is likely the destruction of PAHs within H II regions. Outside of H II regions, $R_{\rm PAH}$ is relatively constant at a value of $\sim 3.43$, indicating either a constant $q_{\rm PAH}$ in the diffuse gas or possibly that we are in a high $q_{\rm PAH}$ regime, where the mapping between $R_{\rm PAH}$ and $q_{\rm PAH}$ begins to flatten (see Section 3.2 for details). Comparing to the Hensley & Draine (2023) models, we can determine what possible values of $q_{\rm PAH}$ correspond to the average value of $R_{\rm PAH}$ we measure in the diffuse gas. A range of possible $q_{\rm PAH}$ values for different single-$U$ radiation field strengths and PAH size distributions are listed in Table 4. We find that our average $R_{\rm PAH,Diffuse}$ of 3.43 corresponds to a $q_{\rm PAH}$ of 5.5% using the Hensley & Draine (2023) dust models assuming the Draine et al. (2021) standard grain distribution, 8.1% assuming the large grain distribution, or 3.91% using the small grain distribution, in a radiation field with $\log_{10}$U=0.0 and $\gamma=0$. We note that fixing $\gamma=0$ (i.e. using single values of $U$ rather than a distribution) will produce the widest range of inferred $q_{\rm PAH}$ values for a given $R_{\rm PAH}$, since there is no contribution to the MIR emission from radiation fields that extends up to high intensities. As can be seen in Draine & Li (2007) Figure 18, contributions from the power-law component can be important at 21 µm. Therefore, Table 4 represents the widest range of $q_{\rm PAH}$ measurements that could yield $R_{\rm PAH}$ $=3.43$. In the remainder of this Section, we will examine any secondary trends in $R_{\rm PAH,Diffuse}$ that we can measure with the full suite of ancillary datasets available.

The fraction of dust in the form of PAHs with $<10^{3}$ C atoms ($q_{\rm PAH}$, Draine & Li, 2007) can be constrained using MIR through FIR observations. From the FIR, the peak of the thermal dust SED provides information on the radiation field intensity. The shape of the FIR SED can also provide constraints on the distribution of radiation fields heating the dust (Dale et al., 2014). With this information on the radiation field, the addition of MIR measurements that sample PAH emission yield a constraint on $q_{\rm PAH}$, since the primarily stochastic emission from PAHs is proportional to the radiation field intensity. The inferred PAH surface density compared with the dust mass surface density constrained by the FIR peak, yields $q_{\rm PAH}$.

The resolution of existing FIR datasets strongly limits the sample of galaxies where $q_{\rm PAH}$ can be resolved at $10-50$ pc resolution to the Local Group. For this reason, it is critical to explore MIR-only indicators of PAH fraction to extend observations beyond the nearest galaxies in order to make progress in understanding the life-cycle of PAHs. In this paper we have explored the MIR-only diagnostic $R_{\rm PAH}$ $=$(F770W_ss+F1130W)/F2100W as an indicator of $q_{\rm PAH}$. We have compared $R_{\rm PAH}$ and $q_{\rm PAH}$ at low resolution for several galaxies with FIR observations that were also covered by PHANGS-JWST and found good correlation. We also explored the relationship between the predicted $R_{\rm PAH}$ to $q_{\rm PAH}$ in the Draine & Li (2007) dust models and their recent updates (Draine et al., 2021; Hensley & Draine, 2023). In general, at dust-mass weighted average radiation field intensities $\lesssim 10^{2}$, which is typical for regions in our galaxies at $10-50$ pc resolution, $R_{\rm PAH}$ is well correlated with $q_{\rm PAH}$, yielding a good diagnostic for changes in the PAH fraction (see also Appendix E). It is worth noting that the mapping between $R_{\rm PAH}$ and $q_{\rm PAH}$ is non-linear and at high values of $q_{\rm PAH}$ $>5$% (super-Milky Way $q_{\rm PAH}$ values), the slope between them becomes flatter, resulting in only small changes in $R_{\rm PAH}$ as $q_{\rm PAH}$ increases. On the other hand, decreases in $R_{\rm PAH}$ are a good diagnostic for revealing locations where the PAH fraction drops, presuming the dust-mass weighted average radiation field is not $U\gg 10^{2}$.

There are several questions about the $R_{\rm PAH}$-$q_{\rm PAH}$ relationship that should be addressed in future studies. First, we have focused on the Draine & Li (2007) model, which makes specific assumptions about the dust grain size distribution. Other dust models make different assumptions about grain size distributions and dust optical properties which can change the interpretation of $R_{\rm PAH}$, particularly if there is a variable amount of “very small grains” (VSGs) that can contribute to the MIR continuum at $\sim$20 µm. Future work with DustEM and the THEMIS model (Compiègne et al., 2011; Jones et al., 2013; Ysard et al., 2024) could explore $R_{\rm PAH}$ as a tracer of small hydrocarbon grain abundance in that model framework, which includes variable VSG abundance.

In addition, in our FIR modeling comparison and exploration of the model predictions for $R_{\rm PAH}$-$q_{\rm PAH}$ mapping, we have assumed a radiation field distribution (specifically, a delta function plus power-law model). This choice directly impacts the interpretation of the MIR emission, since the power-law component (with a distribution of radiation fields extending from $U_{\rm min}$ to $U=10^{7}$ and a slope fixed to $\alpha=2$) contributes MIR continuum emission at $\sim 20$ µm. In models that assume a single radiation field, that MIR emission is often attributed to changes in the grain size distribution (specifically, an increase in very small grains, see Figure III.15 of Galliano, 2022, for a very clear visualization of the differences between these types of models). The applicability of a delta function plus power-law distribution has been explored in a wide variety of contexts for unresolved galaxies or $\sim$kpc scale measurements (Dale et al., 2005; Aniano et al., 2012; Galliano et al., 2018). At a very high spatial resolution, in theory one might expect to reach a scale where a single radiation field may be applicable, but in general on $10-50$ pc scales it seems likely that a distribution of radiation field intensities is present. Little is known about the detailed small scale distribution of radiation field intensities throughout the ISM. Future work investigating this will be critical to improve understanding of the $R_{\rm PAH}$-$q_{\rm PAH}$ mapping and for reconciling differences between models that interpret the MIR continuum as changes in small grain abundance, PAH fraction, and radiation field intensity and distribution.

H II regions crisply stand out in our $R_{\rm PAH}$ maps as regions with low PAH fraction. This observation agrees well with previous studies in MW H II regions (e.g. Cesarsky et al., 1996; Kassis et al., 2006; Povich et al., 2007), with extragalactic observations that have the spatial resolution to isolate H II regions (e.g. Chastenet et al., 2019), and with early results on JWST observations of our sample (Chastenet et al., 2023a; Egorov et al., 2023). Additionally, this follows from the trends between sSFR (measured by 70µm surface brightness) and $q_{\rm PAH}$ measured at lower spatial resolution described in (Aniano et al., 2020). When we separate “nebular” regions using the Groves et al. (2023) catalog, we find that these boundaries encompass nearly all of the low $R_{\rm PAH}$ regions in our maps. We also find that the $R_{\rm PAH}$ weighted by H$\alpha$ shows an even more dramatic decrease in the nebular regions. To confirm this result is valid even if some H II regions were not included in the nebular catalog of Groves et al. (2023), we also examine trends in $R_{\rm PAH}$ as a function of $I_{\rm H{\alpha}}/\Sigma_{\rm H{\sc I}+H_{2}}$, a pixel-based tracer of the amount of ionized gas, and find that $R_{\rm PAH}$ falls dramatically at a value of $I_{\rm H{\alpha}}/\Sigma_{\rm H{\sc I}+H_{2}}$$=37.5$erg s^-1 kpc^-2 (M_⊙pc^-2)^-1, confirming that $R_{\rm PAH}$ decreases in regions dominated by ionized gas.

The good correspondence of the presence of H II regions with the decrease in $R_{\rm PAH}$ across a sample of H II regions with a wide range of ages and luminosities suggests that PAH destruction occurs quickly compared to the H II region lifetime. Indeed, observations of PAH emission in MW H II regions suggest PAH survival times on the order of $\sim$few thousand years (Kassis et al., 2006; Compiègne et al., 2007) inside the ionization front. Recent JWST observations of the Orion Bar photodissociation region (PDR) have emphasized the sharp decrease in PAH emission near the ionization front (Chown et al., 2023; Peeters et al., 2023), corroborating a quick destruction timescale in the photoionized gas.

The exact mechanism of PAH destruction in H II regions is not fully understood. Laboratory studies of PAHs have suggested that fragmentation of PAH molecules is likely to occur when they are irradiated by photons with energies between 8–40 eV, although there is some uncertainty as to what photon energies will primarily fragment PAHs and which energies will lead only to photoionization of the PAHs (Jochims et al., 1996; Zhen et al., 2015). Theoretical calculations suggest that in photoionized gas with H II region-like densities, thermal sputtering would require long timescales ($\sim 10^{7}$ years; Micelotta et al., 2010), which would be difficult to reconcile with the strong correspondence between H II region locations and dips in $R_{\rm PAH}$ across our large sample of regions. Observations of MW H II regions provide evidence for extreme UV (EUV) photons playing an important role in PAH destruction (Povich et al., 2007). Comparison of optical line ratios from MUSE to the regions of low $R_{\rm PAH}$ in the PHANGS-JWST targets also show $R_{\rm PAH}$ is correlated with tracers of the hardness and intensity of the UV field (Egorov et al., 2023). More detailed investigation of the correlation of $R_{\rm PAH}$ and PAH properties judged from spectroscopy with tracers of the ionized gas conditions and radiation field should provide key insights into the destruction pathways.

Because of the resolution of the MUSE and JWST data, our nebular regions always include contributions from both the H II region, surrounding PDR, and depending on the H II region size and galaxy distance and inclination, a substantial contribution from the surrounding “diffuse” gas. This means that the low $R_{\rm PAH}$ signal associated with H II regions is diluted, particularly for small H II regions (see Section 4.1.3 for more discussion). Future work comparing the $R_{\rm PAH}$ maps to high resolution H$\alpha$ from HST observations or Paschen-$\alpha$ from JWST should better resolve the boundary of the H II region enabling less contaminated measurements of the decrease in $R_{\rm PAH}$. We do see, however, that in cases where H II regions are large or the line of sight is dominated by ionized gas (see Figure 10) the decrease in $R_{\rm PAH}$ is more stark. This is likely to be a combination of a better match between the MUSE nebular region and the actual H II region size, and the fact that large H II regions may be big enough to fully puncture through the galaxy disks, giving a line of sight that has little foreground or background contamination from the surrounding neutral ISM. It is worth noting that in the case of $q_{\rm PAH}$ $\sim 0$, $R_{\rm PAH}$ does not asymptote to 0, but something more like $R_{\rm PAH}$ $\sim 0.5$ in the Draine & Li (2007) and Hensley & Draine (2023) models. The nature of the MIR continuum in situations where there are truly no PAHs is not well characterized since few lines of sight with those characteristics have been studied in the tracers needed to constrain the grain population in the MW. Observations with Spitzer and SOFIA of Galactic H II regions NGC 7023 (Croiset et al., 2016) and the Eagle Nebula (Flagey et al., 2011) show little to no PAH emission inside the photodissociation regions (PDRs), but these observations were not the main comparison source for the Draine & Li (2007) and Hensley & Draine (2023) models. This means that at very low $q_{\rm PAH}$, the $q_{\rm PAH}$-$R_{\rm PAH}$ mapping may be more uncertain. Regardless, it is clear that we see the lowest $R_{\rm PAH}$ values in lines of sight that are most dominated by ionized gas. Our observations are consistent with a picture where ionized gas is essentially free of PAHs, as would be expected based on observations of MW H II regions and estimates of the PAH destruction timescale.

The observation that H II regions are holes in the distribution of PAHs in the ISM provides a potential explanation for the observed trends of kpc-scale $q_{\rm PAH}$ and $R_{\rm PAH}$ with SFR and sSFR (e.g. Aniano et al., 2020). As the SFR surface density increases, a larger fraction of the ISM in that kpc is comprised of H II regions, leading to a lower average PAH fraction. Also, because the kpc-scale averages are luminosity-weighted, and the surroundings of H II regions have some of the highest MIR luminosities, low $q_{\rm PAH}$ material will contribute more to setting the kpc-average $R_{\rm PAH}$. An interesting consequence of the SFR or sSFR trends in $R_{\rm PAH}$ may be to exaggerate the observed dependence of $q_{\rm PAH}$ on metallicity. Because low metallicity galaxies tend to have low stellar mass and therefore have higher sSFR (and in particular many of the low metallicity galaxies previously studied in the MIR are starbursts with very high sSFR; Wu et al., 2006; Engelbracht et al., 2008), some part of the observed metallicity trend may in fact be a sSFR trend. Future JWST observations should be able to clearly disentangle these two effects by resolving the $R_{\rm PAH}$ maps of these sources and isolating the H II regions from the surrounding neutral ISM.

One of the strongest observed environmental trends in PAH fraction is its dependence on metallicity (Engelbracht et al., 2005; Madden et al., 2006; Draine et al., 2007; Rémy-Ruyer et al., 2015; Chastenet et al., 2019; Shivaei et al., 2024). Galaxy integrated $q_{\rm PAH}$ measurements suggest a decrease in the PAH fraction at $12+\log$(O/H)$\sim$8.1-8.3 (Draine et al., 2007; Aniano et al., 2020; Galliano et al., 2021). This decline has been attributed to enhanced PAH destruction (Madden et al., 2006; Gordon et al., 2008) or impeded PAH formation (Sandstrom et al., 2010) at low metallicity, among other possibilities. Observations of the Magellanic Clouds have revealed that at 10 pc resolution, where H II regions and diffuse gas can be separated, the diffuse neutral gas PAH fraction shows a sharp decrease between the metallicity of the LMC and SMC, and that the diffuse neutral gas in the LMC has a nearly MW $q_{\rm PAH}$ value of $\sim 5$% (Chastenet et al., 2019).

Understanding the origin of this metallicity trend is critical for using PAH emission to trace SF and for predicting attenuation curves for galaxies across a range of metallicities and redshifts. It also holds key information about the PAH life cycle—a subject that remains poorly understood. As discussed in Aniano et al. (2020), the galaxy-integrated $q_{\rm PAH}$ measurements as a function of metallicity do not provide a clear distinction between a model where the PAH fraction decreases gradually with metallicity or where there is a sharp decrease at some threshold metallicity.

Although our sample does not span a large metallicity range, nor does it reach the $12+\log$(O/H)$\sim$8.1-8.3 metallicities where the PAH fraction becomes substantially lower than MW, our observations of a nearly constant $R_{\rm PAH,Diffuse}$ throughout our sample suggest that metallicity trends are steep, rather than gradual. We see no clear metallicity trends in $R_{\rm PAH,Diffuse}$ over $\sim 0.3$ dex in $12+\log$(O/H) and in fact one of our lowest metallicity galaxies has one of the highest $R_{\rm PAH,Diffuse}$ values. The lack of any systematic metallicity trend may lend some support to the scenario where there is a sharp drop in $q_{\rm PAH}$ at some metallicity. This agrees well with the LMC results (Chastenet et al., 2019) showing a MW-like $q_{\rm PAH}$ in the diffuse neutral gas of that galaxy, despite its low metallicity. An important future direction for understanding PAH fraction as a function of metallicity is to separate $R_{\rm PAH}$ in H II regions and diffuse gas in lower metallicity targets. If the metallicity dependence is sharp, we may expect to see $R_{\rm PAH,Diffuse}$ in galaxies persist at high values until some threshold metallicity.

Outside of the H II regions identified with the nebular catalogs of Groves et al. (2023), we find that $R_{\rm PAH}$ is fairly constant, with an average value of $\sim$3.4. Based on the dust emission models of Hensley & Draine (2023), this corresponds to a $q_{\rm PAH}$ of 5.55%, assuming the “standard” PAH size distribution and a single radiation field with $U=1$. The $q_{\rm PAH}$ calibration in Hensley & Draine (2023), places the MW diffuse ISM at $q_{\rm PAH}$ $=5.9$%, very similar to what we infer in the diffuse ISM of galaxies in our sample. If we instead translate $R_{\rm PAH}$ to $q_{\rm PAH}$ using a distribution of radiation fields with $U_{\rm min}=1$ and $\gamma=0.005$ in the Draine & Li (2007) model, we find that $R_{\rm PAH}$ $=3.4$ translates to $q_{\rm PAH}$ $\sim 4.5-5$ (see Figure 7). Within the Draine & Li (2007) model, the MW diffuse ISM $q_{\rm PAH}$ $=4.55$%, again very similar to what we infer for our galaxies. These comparisons suggest that we are seeing $q_{\rm PAH}$ values in the diffuse neutral gas of our galaxies that is very similar to what is found in the Milky Way. In addition, the consistency of $R_{\rm PAH}$ outside of nebular regions suggests that the diffuse ISM of these galaxies has a fairly uniform $q_{\rm PAH}$. The small degree of environmental variation at $\sim$Solar metallicity in PAH fraction also helps explain observations showing a tight, nearly linear correlation of PAH emission with gas column (Chown et al., 2021; Leroy et al., 2023b). This is similar to what has been found in the Milky Way, where gas column and dust emission are tightly correlated in the high latitude cirrus (Boulanger et al., 1996; Planck Collaboration et al., 2011; Lenz et al., 2017).

The most significant driver of changes to the average value of $R_{\rm PAH,Diffuse}$ across our sample is the specific star formation rate (sSFR), a correlation that we tested both over the full JWST footprint (Figure 9) and within 1 kpc regions (Figure 11). The positive correlation between $R_{\rm PAH,Diffuse}$ and sSFR could be explained if regions with higher sSFR have an overall harder radiation field in the diffuse ISM compared to lower sSFR regions, which would increase $R_{\rm PAH}$ measured in the diffuse gas at a given $q_{\rm PAH}$. Future efforts to model the interstellar radiation field using information from the stellar populations observed by MUSE should provide important insights into these trends. In addition, studies of $R_{\rm PAH}$ in the vicinity of star clusters, where the radiation field can be predicted with stellar population modeling, show trends consistent with $R_{\rm PAH}$ responding to changes in radiation field hardness (Dale et al., 2023b, Dale et al. in prep). Lastly, work with PHANGS-JWST exploring PAH band ratios has shown evidence of the dependence of PAH emission feature strengths on the hardness of the radiation field (Chastenet et al., 2023b; Baron et al., 2024). Accounting for the response of PAH emission to changes in radiation field spectrum will be a critical next step in interpreting $R_{\rm PAH}$ variations.

In addition to the correlation between $R_{\rm PAH,Diffuse}$ and sSFR, we see a slight decreasing trend in $R_{\rm PAH,Diffuse}$ and $f_{\rm{Mol}}$. This decrease could be explained by PAHs coagulating or sticking to larger dust grains in the dense gas. Ysard et al. (2013) suggest that coagulation of dust grains will occur above a density threshold of a few 10³ H cm^-3. Above this threshold, smaller grains may form fluffy aggregates, producing a different dust population than in the lower density diffuse ISM and lowering $q_{\rm PAH}$. This aggregation could potentially cause some of the decrease in $R_{\rm PAH,Diffuse}$ we observe at high $f_{\rm{Mol}}$, where we expect the gas to be the densest. Additionally or alternatively, in regions of dense molecular gas, $A_{V}$ is expected to be highest, blocking UV light that would excite PAHs from reaching these regions. This high attenuation of UV light could also be lowering $R_{\rm PAH}$ even if $q_{\rm PAH}$ remains constant, it softens the radiation field, leading to lower $R_{\rm PAH}$ at a given $q_{\rm PAH}$. Determining whether the main cause of this decreasing trend is coagulation or UV attenuation requires further spectral information not available with this dataset. This trend appears consistent across our full sample with the exception of NGC 1300 and NGC 5068, which show slight increases in $R_{\rm PAH,Diffuse}$ at high $f_{\rm{Mol}}$.

It should also be noted that our isolation of $R_{\rm PAH}$ from diffuse gas is limited by our ability to identify H II regions with MUSE. Some of the decrease observed in $R_{\rm PAH,Diffuse}$ at high $f_{\rm{Mol}}$ could be caused by embedded star-forming regions undetected by MUSE (as found in a small set of cases in Hassani et al., 2023). These embedded regions would decrease $R_{\rm PAH,Diffuse}$ through PAH destruction, as observed in the lower $R_{\rm PAH,Neb}$ measurements above. The magnitude of this effect is likely to be limited by the short duration of this highly embedded phase of star formation (Kim et al., 2023) and the overall excellent correspondence between 21 µm sources and H$\alpha$ (with 92% of F2100W bright sources falling within an identified nebular region Hassani et al., 2023).