The Metallicity Dependence of PAH Emission in Galaxies. I. Insights from Deep Radial Spitzer Spectroscopy

Cory M. Whitcomb; J.-D. T. Smith; Karin Sandstrom; Carl A. Starkey; Grant P. Donnelly; Bruce T. Draine; Evan D. Skillman; Daniel A. Dale; Lee Armus; Brandon S. Hensley; Thomas S.-Y. Lai; Robert C. Kennicutt

doi:10.3847/1538-4357/ad66c8

1. Introduction

The smallest carbonaceous dust grains in the interstellar medium (ISM) emit bright, broad vibrational features in the mid-infrared (MIR). These dust grains are commonly believed to be polycyclic aromatic hydrocarbons (PAHs) composed of dozens to thousands of aromatically bonded carbon atoms (Tielens 2008; Li 2020). The emission from these PAH features alone contributes over 10% of the total infrared (TIR) luminosity of star-forming galaxies (Smith et al. 2007).

Many previous studies across a wide variety of interstellar conditions have found deficits in the fraction of PAH luminosity relative to total dust continuum power at subsolar metallicity, typically interpreted as a marked change in the abundance of dust in the form of PAHs (e.g., Engelbracht et al. 2005; Jackson et al. 2006; Madden et al. 2006; Draine et al. 2007; Smith et al. 2007; Gordon et al. 2008; Muñoz-Mateos et al. 2009; Haynes et al. 2010; Sandstrom et al. 2012; Aniano et al. 2020; Lai et al. 2020; Whitcomb et al. 2020; Zang et al. 2022; Shivaei et al. 2024). Little is known about the origins of this deficit. In the era of JWST, it is crucial that we understand this PAH–metallicity relation (PZR) in order to accurately interpret PAH emission from local highly resolved galaxies to systems at high redshift where metallicities for a given stellar mass may be significantly lower.

A common explanation for the PZR is that smaller PAHs are destroyed by high-energy photons in low-metallicity regions (Madden et al. 2006; Gordon et al. 2008; Egorov et al. 2023). At low metallicity, young stars are hotter and less dust is present to attenuate their hard photons (Massey et al. 2005). The resulting interstellar radiation fields (ISRFs) have a higher ratio of UV-to-optical photons, increasing the stochastic heating rate and peak temperature of PAHs. In general, the heat capacity of dust grains is proportional to their volume, so the smallest PAH molecules may begin to dissociate when they absorb photons of sufficiently high energy (Guhathakurta & Draine 1989; Le Page et al. 2003; Duley 2009; Draine et al. 2021).

In this paper, we investigate the PZR using large, deep, Spitzer spectral maps (5–18 μm) in radial strips across three nearby galaxies—M101, NGC 628, and NGC 2403—with well-studied abundance gradients obtained using auroral line–based measurements (Kennicutt et al. 2003; Berg et al. 2013, 2015; Croxall et al. 2016; Rogers et al. 2021). These MIR spectral mapping data sets are some of the deepest ever observed with the Spitzer InfraRed Spectrograph (Spitzer-IRS) and extend up to ∼15 kpc from galaxy center, enabling variations in the full PAH emission spectrum to be tracked, aside from the weak 3.3 μm band. By studying how the PAH spectrum changes in strength and shape as a function of radius within each galaxy, we can isolate the variations that result primarily from the radial metallicity gradient. For the first time, we model the recovered PZR trends using detailed physical models of PAH emission in the context of varying incident radiation fields and underlying grain size distributions (GSDs).

This paper is organized as follows. In Section 2, we summarize the key physical processes and potential explanations for the observed PZR trends. In Section 3, we describe our data, region definition criteria, and the resulting suite of data products. In Section 4, we describe details of the observed PZR trends. In Section 5, we introduce several physically motivated model scenarios and test them against the observations. In Section 6, we discuss the physical implications of our results and compare them with previous studies. In Section 7, we outline our most significant conclusions.

2. Background and Definitions

As the gas-phase metal abundance of the ISM drops, the dust-to-gas mass ratio decreases toward the background level set by stellar production (Rémy-Ruyer et al. 2014), or lower in galaxies where grain destruction is more efficient than grain growth (Aniano et al. 2020). The resulting lack of dust opacity increases the average hardness 〈h ν〉 and intensity U of the starlight incident on dust grains. Just as important as the heating environment, the distribution of grain sizes plays a key role in the emerging PAH spectrum.

Leading models suggest the distribution of grain sizes results from a complex balance between many opposing forces: injection in carbon-rich asymptotic giant branch (AGB) winds, destruction in supernova shocks and ionized gas, fragmentation of larger grains into smaller units, coagulation and accretion building grain mass in dense clouds, and so on (Dwek & Scalo 1980; Zhukovska et al. 2008; Draine 2009; Micelotta et al. 2010a; Galliano et al. 2018; Choban et al. 2022, 2024; Narayanan et al. 2023). Many of these processes are sensitive to metal abundance. For example, the timescale of gas-phase accretion onto grain surfaces has a strong dependence on the free metal content of the ISM (Zhukovska et al. 2008; Zhukovska & Henning 2013).

One of the earliest Spitzer results on PAHs revealed a striking difference in PAH emission between galaxies with solar metallicity and those below one-quarter solar: over this range, PAH emission declines rapidly (Engelbracht et al. 2005; O'Halloran et al. 2006; Smith et al. 2007; Hunt et al. 2010; Rémy-Ruyer et al. 2015; Shivaei et al. 2017; Chastenet et al. 2019; Aniano et al. 2020). It remains unclear whether the PAH abundance fraction relative to the total dust mass—q_PAH—drops suddenly or smoothly as metallicity decreases. Some studies find q_PAH remains constant at about 4%, then drops suddenly to about 1% after a threshold at about 0.2–0.3 Z_⊙ (Draine et al. 2007; Chastenet et al. 2019; Shim et al. 2023). Others find a continuous decrease in the PAH fraction as a function of metallicity (Galliano et al. 2008, 2018; Rémy-Ruyer et al. 2015; Aniano et al. 2020).

It has been common to study the PZR using integrated galaxy measurements from galaxies across a wide range of metallicity (Engelbracht et al. 2005; Draine et al. 2007; Smith et al. 2007; Engelbracht et al. 2008; Hunt et al. 2010; Shim et al. 2023). This method introduces complications because the radiation environments and star formation histories vary between galaxies. Resolved studies within individual galaxies have also found similar deficits in relative PAH-to-dust luminosity as metallicity decreases (Lebouteiller et al. 2011; Khramtsova et al. 2014; Chastenet et al. 2023). Within galaxies, the main properties that vary with radius on scales larger than 100 pc are the stellar density, stellar age, and gas-phase metal abundance. As the density of stars decreases at larger radii, the average intensity U of the ISRF decreases as well. The average age of stars and the ISM density also decrease at larger radii (García-Benito et al. 2017; Dale et al. 2020), increasing the average hardness 〈h ν〉 of the ISRF.

The effect each of these properties of the radiation environment has on the PAH spectrum can be explicitly modeled with the new suite of theoretical spectra from Draine et al. (2021, hereafter D21). The D21 models are based on the simplifying assumption that PAHs can be idealized as nearly spherical grains of pure graphite characterized by a single size parameter and a binary ionization state. This is sufficient to produce MIR emission features that match observations of galaxies. Other models use spectra of specific PAH molecules from quantum-chemical calculations; however, this is currently only feasible for a small fraction of the PAH size distribution (Rigopoulou et al. 2021). Because we are interested in exploring the spectrum produced by the full range of PAH sizes exposed to various radiation environments, the D21 models provide an excellent framework for comparison to the observations.

Resolved PZR studies probing the centers of individual star-forming spiral galaxies using photometric tracers of PAH emission have been conducted with JWST imaging. These have well-constrained metallicity gradients in the radial direction, and the fraction q_PAH is found to decrease relatively smoothly (Chastenet et al. 2023). Most previous studies on the PZR only make use of photometric PAH tracers, primarily at 8 μm to study the strongest PAH feature at 7.7 μm. This feature is primarily emitted by ionized PAHs with ∼100 carbon atoms (N_C), resulting in the common explanation that the smallest PAHs are being preferentially destroyed in low-metallicity environments. The full PAH spectrum has contributions from both ionized and neutral grains with ∼20–10⁶ N_C. A full accounting of the properties of the small grain population requires spectroscopic access to all the main PAH emission bands.

In addition to changes in the overall PAH luminosity, we can track several key variations in the PAH spectrum itself as metallicity decreases. Based on the theoretical work of Li & Draine (2001), Draine & Li (2007), and D21, each PAH feature traces a specific regime of the grain size and ionization distributions of the local PAH population. In general, small PAHs emit more efficiently at shorter wavelengths and large PAHs emit more efficiently at longer wavelengths. Specifically, the approximate peaks are as follows (D21; see Figure 6): the 17 μm feature is emitted most efficiently by PAHs with ∼2000 N_C with no ionization preference, 11.3 μm by neutral PAHs with ∼400 N_C, 7.7 μm by ionized PAHs with ∼100 N_C, and 6.2 μm by ionized PAHs with ∼70 N_C. The final major PAH feature at 3.3 μm is emitted at peak efficiency by neutral PAHs just above the minimum survivable PAH size adopted in D21 of 27 N_C. Since this work uses Spitzer spectroscopy from 5 to 38 μm, we do not have 3.3 μm coverage. Instead, we use the best-fit models to predict the metallicity trends in this final PAH feature.

Smith et al. (2007) found that the ratio of the 17–11.3 μm features drops significantly as metallicity declines. This was attributed to enhanced formation of grains larger than ∼400 N_C at higher metallicities. Sandstrom et al. (2012) found a similar result: the 17 μm feature is significantly weaker, but still detected, in low-metallicity regions of the Small Magellanic Cloud (SMC) compared to the higher-metallicity galaxies of the Spitzer Infrared Nearby Galaxies Survey (SINGS; Smith et al. 2007). They found the 7.7 μm feature is also significantly weaker and the 6.2 and 11.3 μm features carry an increased fraction of the PAH power. This suggests PAHs in the low-metallicity SMC are smaller and more neutral. Hunt et al. (2010) found the fractional PAH power in the 8.6 μm and 11.3 μm features is increased for low-metallicity blue compact dwarf galaxies compared to the SINGS galaxies, while the 6.2 μm and 7.7 μm features remain constant. They interpret this as evidence for a decrease in the PAH population smaller than 100 N_C in low-metallicity galaxies compared to high-metallicity systems.

Studies of the PZR in galaxies have resulted in many different interpretations, and several questions remain unanswered. Some attribute the PAH deficit to timing: a galaxy with low metallicity and a young stellar population simply has not had enough time to build up a significant reservoir of PAHs (Galliano et al. 2008). However, most agree that the relative PAH-to-dust luminosity drops as a result of (or correlated with) the decreased gas-phase metal abundance (O'Halloran et al. 2006; Draine et al. 2007; Gordon et al. 2008; Hunt et al. 2010; Shivaei et al. 2017; Aniano et al. 2020; Shivaei et al. 2024).

Interpretation of the observed PZR trends is further complicated by the many unknowns about how small grains are produced and survive in the ISM. It is still unclear whether the trends are a result of changes to the initial PAH GSD from stellar injection or changes in formation and processing pathways in the ISM. What the dominant sources of grains that emit as PAHs are and how their formation mechanisms depend on metal content are areas of ongoing research (Zhukovska et al. 2008; Galliano et al. 2018). It has been theorized that the smallest dust grains should be more susceptible to destruction by ionizing radiation in low-metallicity environments due to their reduced heat capacities (Micelotta et al. 2010a, 2010b; D21). A deficiency of PAH emission has been observed in extreme environments such as H ii (Cesarsky et al. 1996; Kassis et al. 2006; Povich et al. 2007; Gordon et al. 2008; Chastenet et al. 2019; Egorov et al. 2023), likely due to strong radiative processing. Therefore, we expect in some cases that it is possible for PAHs to be destroyed while larger dust grains survive. It is likely that grain formation and destruction processes shape the observed spectrum, and these may be more or less dominant in different environments and galaxies.

3. Data

The primary data products used in this study are cospatial Spitzer-IRS Short-Low (SL1 7.5–14.7 μm, SL2 5.2–7.6 μm, and SL3 7.4–8.7 μm) and Long-Low (LL1 20.6–38.4 μm and LL2 14.2–21.1 μm) spectroscopic maps. We extracted spectra from these maps in ∼0.43–1.7 kpc apertures across three nearby galaxies from the PAH-BIGMAP program 20518: M101 (NGC 5457), NGC 628, and NGC 2403 (see Table 1). These maps extend radially outward from the center of each galaxy at an orientation designed to avoid bright H ii regions. We also include eight Spitzer spectral observations of H ii regions in M101 from Gordon et al. (2008), rereduced in a consistent manner.

Figure 1 shows the position of each spectroscopic map on a three-color image of each galaxy. These images are composed of Infrared Array Camera (IRAC) 3.6 μm in blue, IRAC 8 μm in green, and Multiband Imaging Photometer for Spitzer (MIPS) 24 μm in red (see Section 3.2). Note the shift from green to red with increasing radius in each galaxy. This suggests the PAH emission at 8 μm is decreasing faster than the dust continuum emission at 24 μm.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Three-color images for M101 (top), NGC 628 (left), and NGC 2403 (right). Red is MIPS 24 μm, green is IRAC 8 μm, and blue is IRAC 3.6 μm. All images are convolved to match PACS 160 μm resolution, except those for NGC 2403, which are convolved to a 25'' Gaussian. The footprint of the radial spectroscopic maps studied in this work is shown in cyan, and the eight additional apertures for the H ii region maps in M101 are shown as black squares.
Download figure:
Standard image High-resolution image

3.1. Spectroscopic Data Reduction

All spectral maps were reduced with CUBISM (Smith et al. 2007), using background frames created from both the dedicated background spectra associated with each observation set and emission-free regions on the outskirts of the radial strip maps. Automatic and by-eye pixel selection was employed to remove bad pixels and mitigate striping in the final cube.

Slight mismatches in the surface brightness of aperture-matched spectra were found in faint regions for suborders SL1, SL2, and SL3. These were determined to arise at low surface brightness levels from light contamination, which varies spatially on the detector and over time. We corrected the effect with a row-by-row boxcar average of interorder pixels across time-adjacent background-subtracted detector images. This process is implemented with the program sl_io_correct (Sandstrom et al. 2012). Correction using interorder light results in a significant improvement in the flux agreement between SL1, SL2, and SL3 intensities at the wavelengths where they overlap for low surface brightness spectra.

The suborders of SL and LL were stitched by averaging the brightness at wavelengths where they overlap. The combined LL spectrum was then scaled to match the SL spectrum. This was done with a multiplicative factor or an additive offset based on the brightness of the spectrum. After removing backgrounds, residual baseline variations in dark regions of the cube remained below 1 MJy sr⁻¹. For faint spectra (i.e., average brightness <0.75 MJy sr⁻¹), we employed an additive offset to bring spectral segments into agreement. For brighter spectra, we employ a multiplicative factor. We find this process results in good agreement between SL and LL for all spectra in our sample and all offsets and factors are small, typically ∼0.1 MJy sr⁻¹ (additive) and ∼1.1 (multiplicative).

3.2. Infrared Photometric Data

We also compiled infrared photometric data from 8 μm through 160 μm in order to derive the TIR surface brightness in each of our apertures. The 8 μm data from the IRAC, the 24 μm data from the MIPS, and the 70, 100, and 160 μm data from the Photodetector Array Camera and Spectrometer (PACS) of the Herschel Space Observatory were retrieved from the Key Insights on Nearby Galaxies: A Far-Infrared Survey with Herschel (KINGFISH; Kennicutt et al. 2011) and NGC 2403 from the Very Nearby Galaxies Survey (Bendo et al. 2012). NGC 2403 was never observed with PACS 100 μm, so data at 250 μm from the Spectral and Photometric Imaging Receiver (SPIRE) of Herschel were included so the spectral energy distribution (SED) could be modeled. All photometric data were convolved to the point-spread function (PSF) of the lowest-resolution images used for SED determination—PACS 160 μm, FWHM 11 farcs 2, and a Gaussian PSF with FWHM 25'' for NGC 2403.

3.3. Region Definition

The spectral cubes for M101 and NGC 628 were convolved to the resolution of PACS 160 μm (11 farcs 2) with kernels produced by the STINYTIM program (Sandstrom et al. 2009). Since NGC 2403 has no PACS 100 μm data, it was necessary to convolve the cubes to a Gaussian PSF of FWHM 25'' in order to include longer-wavelength far-infrared photometry in the SED fit. The convolution process (when the image is padded with zeros) also dims the edges of each spectral cube where the convolution kernel falls off the coverage footprint. We trimmed the edges of our spectral cube by 5'' on all sides to preserve as much area as possible while ensuring appropriate data coverage. This 5'' is approximately the size of one pixel in the LL cubes and just less than two pixels in the SL cubes.

Extraction apertures were defined to span the entire overlap region between the SL and LL footprints. We tiled each galaxy with rectangular apertures with a minimum side length of the limiting PSF of the supplemental photometric data (11 farcs 2, 25'' for NGC 2403). An iterative algorithm was used to define apertures based on the continuum brightness at MIPS 24 μm to ensure sufficiently high signal-to-noise PAH detections. The final apertures are shown in Figure 2. The initial apertures for this process were square and spanned the width of each radial strip. If the average surface brightness of the aperture was above 0.25 MJy sr⁻¹ at MIPS 24 μm, then the square was split into a vertical pair and a horizontal pair. Next, we checked whether the vertical or horizontal pair was brighter. If the brightest pair was at least 10% brighter than the original square, or above 1 MJy sr⁻¹, then the algorithm continued on the pair. The pair was then broken into two squares, and the same process was repeated until the minimum size was reached.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Three-color images from slices of the Spitzer-IRS SL1 spectral cubes for M101 (top), NGC 628 (left), and NGC 2403 (right). Red is the surface brightness at the PAH 12.6 μm peak, green is at the PAH 11.3 μm peak, and blue is at the PAH 7.7 μm peak. The apertures defined in Section 3.3 are shown in cyan.
Download figure:
Standard image High-resolution image

From all radial strips there are a total of 160 regions that meet the minimum continuum brightness threshold: 85 in M101, 64 in NGC 628, and 11 in NGC 2403. Including the eight H ii regions, a total of 168 regions with well-measured spectra were analyzed.

3.4. MIR Spectral Fitting

The final extracted Spitzer spectra (5–38 μm) were fit with the IDL program PAHFIT (Smith et al. 2007). From PAHFIT we obtained the integrated intensities of all emission features and their propagated uncertainties. We do not expect significant attenuation in the spectra of the spiral galaxies in our sample, so attenuation was not included in the fit. The subfeatures of major PAH complexes were combined with PAHFIT to produce combined power in the 7.7, 11.3, 12.6, and 17 μm bands. We further combined all PAH features (5–18 μm) into ΣPAH using PAHFIT.

The H ii region spectral cubes are smaller and sparsely sampled by the slit, so we did not convolve them to a common resolution. Instead, we scaled our derived PAH luminosities in these apertures by the ratio between their average native IRAC 8 μm photometry and their average IRAC 8 μm when convolved to PACS 160 μm resolution. The median ratio is 1.94, the mean is 1.86, and the standard deviation is 0.25, so the PAH luminosity of each H ii region was scaled down by about a factor of two, with no change in the band ratios.

3.5. Derived Quantities

From the photometric data we derived TIR luminosity (3–1100 μm, TIR) for each region in M101 and NGC 628 using the calibration presented in Galametz et al. (2013):

$\begin{eqnarray}\begin{array}{rcl}\mathrm{TIR} & = & 2.051\times \nu {I}_{\nu }(24\mu {\rm{m}})+0.521\times \nu {I}_{\nu }(70\mu {\rm{m}})\,\\ & & +0.294\times \nu {I}_{\nu }(100\mu {\rm{m}})+0.934\times \nu {I}_{\nu }(160\mu {\rm{m}}),\end{array}\end{eqnarray} \tag{ 1 }$

where the constants were taken from their Table 3 and I_ν (λ) is the surface brightness at λ microns in megajanskys per steradian.

There are no data at 100 μm for NGC 2403, so we used the calibration from Galametz et al. (2013) that does not depend on this band:

$\begin{eqnarray}\begin{array}{rcl}\mathrm{TIR} & = & 2.126\times \nu {I}_{\nu }(24\mu {\rm{m}})\,+0.670\times \nu {I}_{\nu }(70\mu {\rm{m}})\\ & & +1.134\times \nu {I}_{\nu }(160\mu {\rm{m}}),\end{array}\end{eqnarray} \tag{ 2 }$

where the notation is the same as in Equation (1).

Table 1. Spitzer-IRS Observations

Target Name	R.A. (J2000)	Decl. (J2000)	AORID
M101
SLmap_01	14^h03^m0190	+54^d17^m010	14800128
SLmap_00	14^h03^m1050	+54^d20^m108	14799872
SLmap_02	14^h02^m5330	+54^d13^m512	14800384
SL BG	14^h03^m2650	+54^d37^m350	14800896
LLmap	14^h03^m0190	+54^d17^m010	14799616
LL BG	14^h03^m2650	+54^d37^m350	14800640
NGC 628
SLmap_00	1^h36^m4810	+15^d43^m562	14802432
SLmap_01	1^h36^m4330	+15^d46^m144	14802688
SL BG	1^h36^m5413	+15^d37^m590	14803200
LLmap	1^h36^m4570	+15^d45^m053	14802944
LL BG	1^h36^m5413	+15^d37^m590	14803456
NGC 2403
SLmap_01	7^h36^m4470	+65^d35^m237	14801408
SLmap_00	7^h36^m2890	+65^d33^m245	14801152
SL BG	7^h36^m2969	+65^d27^m276	14802176
LLmap	7^h36^m3680	+65^d34^m241	14801664
LL BG	7^h36^m2969	+65^d27^m276	14801920

Download table as: ASCII Typeset image

We determined metallicities (12 + [O/H]) for our extracted regions in M101, NGC 628, and NGC 2403 using the radial gradient equations presented in Rogers et al. (2021). We defined solar metallicity as Z_⊙ ≡ 12 + [O/H] = 8.69 (Asplund et al. 2009), and we assumed that O/H is a good tracer of the abundance of all metal species. Skillman et al. (2020) found the relative C/O abundance decreases radially in M101 by about 0.4 dex over two effective radii (R_e, see Table 2). Therefore there is likely even less carbon than expected from a linear scaling with the oxygen abundance gradient of each galaxy. The parameters of the metallicity gradient equation for each galaxy are shown in Table 2, where the metallicity is largest at the center (Z_o) and decreases with a constant slope (m) as a function of deprojected galactocentric radius (R): Z = Z_o − m × (R/R_e).

Table 2. Galaxy Information and Metallicity Gradients

Property	M101	NGC 628	NGC 2403
Central abundance (12 + [O/H])	8.70 ± 0.04	8.65 ± 0.05	8.55 ± 0.04
Slope m (dex/R_e)	0.17 ± 0.02	0.09 ± 0.03	0.09 ± 0.03
R_e (arcsec)	197.6	95.4	178.0
Distance (Mpc)	7.4	7.2	3.2
Inclination (deg)	18	5	63
Position angle (deg)	39	12	124

Note. Values adopted from Berg et al. (2020) and Rogers et al. (2021).

Download table as: ASCII Typeset image

We assumed metallicity is purely radial with no small-scale or azimuthal variations. Kreckel et al. (2019) found the metallicity of nearby spiral galaxies varies by up to 0.05 dex on scales of 120 pc, but the radial dependence dominates. This is comparable to the scatter in the metallicity gradient fits employed here. Including a 0.05 dex uncertainty for all metallicity values would not significantly alter the results of this work.

3.5.1. Dust and Galaxy Property Maps

The properties of dust and local environmental conditions can be derived from infrared photometric data from IRAC 3.6 μm, IRAC 4.5 μm, IRAC 5.8 μm, IRAC 8 μm, MIPS 24 μm, PACS 70 μm, PACS 100 μm, and PACS 160 μm. This was done by fitting the SED with the Draine & Li (2007) model using the methodology of Aniano et al. (2020). NGC 2403 was not observed with PACS 100 μm, so an additional far-infrared band, SPIRE 250 μm, was necessary to accurately fit the SED model. The SED models produce resolved maps of quantities such as the fraction of dust mass in the form of PAHs with N_C < 10³ (q_PAH) and the mass-weighted average starlight energy density per unit frequency that is illuminating the local dust ( $\bar{U}$ ). We focus primarily on the maps of these two quantities in this analysis. Values for each region are extracted from these maps using the matched apertures defined in Section 3.3.

4. Observational Results

Figure 3 shows three example spectra from the M101 radial strip. Comparing them, we observe that the 17 μm feature and the brightest feature at 7.7 μm are significantly reduced at lower metallicity. The remaining PAH power is in the 11 μm complexes and the 6.2 μm feature, which is brighter at 0.4 Z_⊙ than at Z_⊙.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Example IRS spectra from M101. The spectrum of the center (red), threshold metallicity (green), and edge regions (blue) from the M101 radial strip. Each spectrum is normalized by the surface brightness at 11.3 μm. The PAHFIT result for each spectrum is shown in gray along with each component of the fit: stellar continuum (cyan), dust continuum (orange), line emission (black), and PAH emission (magenta). Note the relative shift in power from long-wavelength PAH features to short as metallicity drops from Z_⊙ to 0.4 Z_⊙.
Download figure:
Standard image High-resolution image

4.1. PAH to Total Dust Emission Metallicity Relation

The ratio of all PAH emission (ΣPAH) to the TIR brightness is shown in the top panel of Figure 4. The ratio ΣPAH/TIR rises from about 10% to 16% as metallicity decreases from Z_⊙ to $\sim 2/3$ Z_⊙ and then drops below 2% by $\sim 1/4$ Z_⊙. We refer to this threshold metallicity where the drop in ΣPAH/TIR occurs as Z_th ≡ 0.625 Z_⊙ (12 + [O/H] = 8.5). The best-fit relation between ΣPAH/TIR and metallicity has the form

$\begin{eqnarray}\displaystyle \frac{{\rm{\Sigma }}\mathrm{PAH}}{\mathrm{TIR}}( \% )=\left\{\begin{array}{ll}15.1\,\pm \,2.0 & \,\,Z\gt {Z}_{{\rm{th}}}\\ (15.1\pm \,0.2){\left(\displaystyle \frac{Z}{\,{Z}_{{\rm{th}}}}\right)}^{2.6\,\pm \,0.1} & \,\,Z\leqslant {Z}_{{\rm{th}}}\end{array}\right.,\end{eqnarray} \tag{ 3 }$

where the fit at metallicities above Z_th is the uncertainty-weighted median with bulge light correction applied (see Section 4.1.1) and below Z_th the fit is a power law. The value of Z_th is determined by requiring the fit to be continuous at Z_th to within 1σ. The H ii regions in the mid-disk are notably offset from the overall trend; their ΣPAH/TIR is somewhat lower than expected based on their metallicity.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Top: PAH emission deficit at low metallicity. Ratio of all PAH emission to all other dust emission (ΣPAH/TIR) as a function of metallicity. The hexagons and their error bars indicate the mean and standard deviation in six bins of width 0.1 Z_⊙. The gray bin includes only the points below 0.4 Z_⊙. The pink bin falls within the bulge radius of M101 and indicates the mean after a correction for reddened ISRF is applied (see Section 4.1.1). The solid black lines indicate the median above Z_th and a power-law fit below. Bottom: correlation between ΣPAH/TIR and modeled PAH mass abundance fraction (q_PAH). The dashed lines indicate the expected trend from D21 models for three incident ISRFs with U = 1: M31 bulge (red), mMMP (black), and a low-metallicity starburst (cyan). The dotted cyan line indicates the low-metallicity spectrum with U = 100. The solid black line indicates a linear fit with intercept fixed at zero. The hexagons indicate the mean and standard deviation in bins of width 0.5% in q_PAH. The gray bin includes only the points below q_PAH = 1%.
Download figure:
Standard image High-resolution image

The bottom panel of Figure 4 shows the strong correlation between ΣPAH/TIR and q_PAH, where the linear fit has the form

$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Sigma }}{\rm{PAH}}}{\mathrm{TIR}}=(3.69\,\pm \,0.02)\times {{q}}_{{\rm{PAH}}}.\end{eqnarray} \tag{ 4 }$

Also indicated in this figure is the expected ΣPAH/TIR for a given q_PAH from the D21 models. We use their "standard" GSD and rescale it such that the resulting integral for N_C < 10³ gives a desired q_PAH. A spectrum is then produced from each GSD, and ΣPAH/TIR is calculated by PAHFIT and Equation (1). We leave U fixed at 1 and repeat this using three different ISRFs to plot the three dashed lines: the 10 Myr Z_⊙/20 SED from Eldridge et al. (2017), the solar neighborhood SED from Mathis et al. (1983) as modified by D21 (hereafter the mMMP SED), and the M31 bulge SED from Groves et al. (2012). An additional dotted line is also shown where we have increased U to 100 and used the 10 Myr Z_⊙/20 SED. These lines indicate the general effect that increasing the ISRF hardness and intensity has on the PAH and dust luminosities.

4.1.1. Bulge ISRF Correction

Following the work of Draine et al. (2014), we can interpret the decrease in ΣPAH/TIR as metallicity increases from 0.9 Z_⊙ to Z_⊙ as a result of the increasing bulge fraction at low galactocentric radii in M101. NGC 628 and NGC 2403 do not have significant bulges (Fisher & Drory 2010). PAHs in the bulge are excited by an ISRF dominated by older stars, so they emit less compared to PAHs in the disk where the ISRF is more dominated by UV and optical photons from young stars. This effect is less significant for larger dust grains, so ΣPAH/TIR drops as the bulge fraction increases. The red dashed line in Figure 4 shows this using the D21 models: a given q_PAH illuminated by the M31 bulge SED results in ∼50% less ΣPAH/TIR than the same q_PAH illuminated by the solar neighborhood mMMP SED.

Figure 5 shows the $\bar{U}$ values extracted from our apertures. Based on this, the disk of M101 has $\bar{U}$ ≈ 1, and it begins to increase due to the bulge component of the galaxy closer to the center. We consider the bulge radius to be where $\bar{U}$ begins to increase above 1 (about 60'' from galaxy center, approximately where metallicity is 0.9 Z_⊙). Within this radius is also where the drop in ΣPAH/TIR and q_PAH occurs for regions in M101.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** $\bar{U}$ as a function of galactocentric radius. Values detected below 3σ are set to zero.
Download figure:
Standard image High-resolution image

We used the equations from Draine et al. (2014) to correct q_PAH by accounting for the softer radiation field in the bulge of M101. Again using the $\bar{U}$ maps we calculated U_bulge for all regions within the bulge radius as U_bulge = $\bar{U}$ − U_disk with U_disk = 1. However, we find that this shifts points within the bulge radius from 1% below the median q_PAH between Z_th and 0.9 Z_⊙ (4%) to 1% above the median. A correction that brings these points to share the median can be obtained by a small modification to Equation (21) from Draine et al. (2014). The parameter A is tied to the hardness of the radiation field, with higher values corresponding to harder radiation. For instance, the A value for the bulge of M31 is 1.95 and the value for the disk is 4.72. We find changing only the value of A_bulge from 1.95 to 3.85 brings the points in the bulge radius to agree with the median of the points outside this radius (and less than Z_th) in both ΣPAH/TIR and q_PAH. This brings A_bulge closer to A_disk = 4.72. This change is reasonable since M101 has a weak, actively star-forming bulge, while the bulge of M31 considered in Draine et al. (2014) has over 3 times the metal abundance, significantly more dust attenuation, and an older stellar population dominating the ISRF.

4.2. PAH Band Ratios

To better diagnose what is happening with the PAHs, we plot short-to-long wavelength band ratios as a function of metallicity in Figure 6. The top panel shows the ratio PAH 11.3/17 μm is constant until Z_th; then it rises steeply as metallicity drops. Since the 17 μm feature originates from the largest neutral and ionized PAHs, this could be evidence that the largest PAHs (> 2000 N_C) are preferentially depleted at low metallicity. This rise could also be a result of increasing ISRF hardness as metallicity decreases below Z_th, causing a shift in PAH luminosity from longer to shorter wavelength bands (Draine et al. 2021).

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** As a function of metallicity: (top) PAH 11.3/17 μm ratio, (second) PAH 7.7/11.3 μm ratio, (third) PAH 6.2/7.7 μm ratio, (bottom) PAH 6.2/17 μm ratio. Hexagons indicate the mean and standard deviation in six bins with width 0.1 Z_⊙. The vertical line indicates the threshold metallicity Z_th. Ratios of features that were detected below 3σ are not included.
Download figure:
Standard image High-resolution image

The second panel shows the ratio PAH 7.7/11.3 μm decreases away from the threshold metallicity on either side. The 7.7 μm feature is emitted most by ionized PAHs with ∼100 N_C, and the 11.3 μm feature is from neutral PAHs with ∼400 N_C. As metallicity decreases from Z_⊙ to Z_th the PAHs may be getting smaller or more ionized; then below Z_th they may become larger or less ionized. It is also possible that the harder ISRF at lower metallicity is raising the fractional PAH luminosity at short wavelengths and resulting in an increase in the PAH 7.7/11.3 μm ratio; however, this could not explain why the ratio then decreases below Z_th.

The third panel shows the ratio PAH 6.2/7.7 μm is mostly constant from Z_⊙ to Z_th and then increases at lower metallicities. The increase is steeper for the M101 strip regions than for those from NGC 2403 and the M101 H ii regions. The 6.2 μm feature is from ionized PAHS with ∼70 N_C so the increase at low metallicity could be a result of PAHS becoming smaller below Z _th . It is also possible that the increasing ISRF hardness is responsible for the increase in the PAH 6.2/7.7 μm ratio.

The final panel shows the ratio PAH 6.2/17 μm rises steadily as metallicity decreases. This again indicates that the PAHs are becoming smaller on average. Alternatively, the PAHs may remain the same and the ISRF hardness is increasing, causing the PAHs to emit more at shorter wavelengths.

Overall, the trends in these four PAH band ratios suggest the following. PAH 7.7/11.3 μm and PAH 6.2/17 μm rise between Z_⊙ and Z_th, while PAH 6.2/7.7 μm and PAH 11.3/17 μm are constant. This suggests the effect of increasing ISRF hardness is not significant at these metallicities, or else all ratios would show a steady increase. The observed increase in PAH 6.2/17 μm by about an order of magnitude could be interpreted as an increase in the amount of PAHs with ∼70 N_C relative to PAHs with ∼2,000 N_C. Similarly, the observed increase in PAH 7.7/11.3 μm by about 40% could be interpreted as an increase in the amount of PAHs with ∼100 N_C relative to PAHs with ∼400 N_C.

It is also possible that the increase in PAH 7.7/11.3 μm is evidence of PAHs becoming more ionized. This is the only ratio shown that is strongly dependent on the ionization distribution of the PAHs. The ionization fraction of PAHs depends on the conditions of the gas and the grain size, with smaller PAHs being more neutral on average than larger PAHs (D21). If PAHs are becoming smaller then the overall fraction of neutral PAHs will increase as well.

Below Z_th, again all ratios except for PAH 7.7/11.3 μm show the same trend of increasing with decreasing metallicity. The three ratios that increase at low metallicity suggest PAHs overall are becoming smaller and/or hotter. If changes in the radiation field hardness were the primary driver, we would expect the PAH 7.7/11.3 μm ratio to also increase below Z_th. This is not what is observed, suggesting instead that there are also changes in the PAH size distribution and overall ionization fraction.

The bottom panel of Figure 7 shows the differential change in each PAH feature relative to the total PAH luminosity compared to the value at solar metallicity. This figure shows the shifts in power among the PAH features. As metallicity drops from Z_⊙ to Z_th, the PAH power shifts out of the longer wavelength features and into the short wavelength features. When metallicity drops below about Z_th the power begins to shift out of the 7.7 μm PAH feature and toward the 11–13 μm complexes. The most dramatic changes are in the 6.2 μm and 17 μm feature. The 6.2 μm feature increases about 60% from Z_⊙ to 1/2 Z_⊙. The 17 μm feature decreases about 75% over this same metallicity range.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Top: the observed metallicity effect on PAHs from the ratio of each major PAH feature to total PAH luminosity. The gray vertical line indicates the threshold metallicity Z_th. Bottom: fractional change in PAH/ΣPAH relative to the value at Z_⊙. Relative to the total PAH luminosity, the power shifts to shorter wavelengths; then after Z_th, the 7.7 and 17 μm features decrease, the 11.3 μm feature is unchanged, and the 6.2 μm feature increases.
Download figure:
Standard image High-resolution image

5. Modeled PAH and Dust SED Analysis

To interpret the changes in band ratios, we created models of PAHs with various size distributions and illuminated by different radiation fields using the theoretical PAH spectra from D21. By comparing the variations in q_PAH, ΣPAH/TIR, and PAH/ΣPAH, we investigated what changes in PAH size and ISRF can reproduce our observed trends. D21 produced model PAH spectra based primarily on four input parameters: incident radiation field spectrum, incident radiation intensity (U), PAH GSD, and a binary PAH ionization distribution (either neutral or cationized). The standard GSD is described in D21 as a double-peaked log-normal distribution with peak locations characterized by grain sizes a₀₁ = 4 Å and a₀₂ = 30 Å (see Figure 8). Based on the agreement between the linear fit in the bottom panel of Figure 4 and the D21 model using the mMMP ISRF with U = 1, we assumed these parameters are constant in each of the scenarios described below other than the ISRF hardness scenario. For simplicity, we used the standard ionization distribution as a function of grain size. When we generated an infrared dust SED from the D21 models, we also included the emission from larger "astrodust" grains in order to get an accurate measure of TIR.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Size distributions of PAH feature emitting grains. Shown are the size distributions as described in D21 with several values of a₀₁. Note that changing a₀₁ shifts the location of the first peak but the integral remains constant. The photodestruction parameter a_min is shown in red; below this value all PAHs are deleted from the distribution without preserving the integral.
Download figure:
Standard image High-resolution image

For each spectrum we generated from the D21 models, we measured all the PAH features using PAHFIT (Smith et al. 2007) to be consistent with the PAH feature strengths measured in the observed spectra. We used an updated version of the PAHFIT model that includes the 3.3 μm PAH feature (Lai et al. 2020). From the PAHFIT results we derived the sum of all PAH emission (ΣPAH) as well as the ratio of every PAH band relative to ΣPAH (PAH/ΣPAH). TIR was determined by integrating the generated spectrum from 3 to 1100 μm. The value of q_PAH for every spectrum is calculated by the ratio of PAH mass—from integration of the input GSD for grains with N_C < 10³ (∼13.4 Å)—relative to the total astrodust mass.

5.1. Model Scenarios

We considered three realistic scenarios that could reproduce the major observed PZR trends: q_PAH, ΣPAH/TIR, and the main PAH/ΣPAH ratios, all as a function of metallicity. The "ISRF hardness" model fixes the size distribution and varies the hardness of the incident ISRF as metallicity decreases. The "photodestruction" model simulates bottom-up destruction of PAHs by increasing the minimum grain size as metallicity decreases. Finally, the "inhibited growth" model mimics a reduced accretion rate of carbon from the gas phase by decreasing the average PAH size and the overall abundance of carbonaceous grains.

5.1.1. ISRF Hardness Model

First, we test whether changes in the incident ISRF SED alone can explain the observed PZR trends. For this test we use the two ISRF hardness extremes: the bulge SED of M31 and a 10 Myr Z_⊙/20 SED (low-Z). At high metallicity (at galaxy center), we assume the ISRF is weighted as 1:0 M31:low-Z, and at the lowest metallicity (at galaxy edge) we assume the weighting is 0:1 M31:low-Z with a linear interpolation between these extremes. The low-Z SED is characteristic of significantly lower metallicity than we expect at the edges of our three galaxies, and the M31 bulge SED is likewise significantly softer than we expect in the centers of our galaxies. Therefore we expect that the resulting changes to the PAH ratios from this ISRF hardness model represent an upper bound and more realistically the radial change in radiation hardness for our three galaxies will be smaller than suggested by this model.

Figure 9 shows the effect of changing the ISRF hardness on the ratios PAH/ΣPAH: trends broadly similar to the observations are produced, but with a magnitude that is insufficient by a factor of 3–5. Since the ratios from this model represent an upper bound, this indicates that the increased radiation hardness at low metallicity may play a role in the observed PZR trends, but the effect is not strong enough to be the primary driver. This model is also incapable of producing the observed drop in q_PAH and ΣPAH/TIR as metallicity decreases below Z_th. When only the incident radiation field hardness is increased, PAHs of all sizes emit more, so ΣPAH increases. Due to the difference in heat capacities between PAHs and larger dust grains that emit mid- and far-IR continuum, ΣPAH increases faster than TIR in harder radiation environments. Since q_PAH is derived from the GSD and we do not change it in this scenario, q_PAH remains constant regardless of the incident ISRF. In practice q_PAH is often inferred from the 8 μm brightness, and, similar to ΣPAH/TIR, this increases faster than the total dust luminosity as the incident radiation field becomes harder. Therefore changing the ISRF hardness alone cannot explain all observed PZR trends.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Fractional PAH power as a function of metallicity for the ISRF hardness model (Section 5.1.1). The binned data from Figure 7 are shown. Also shown is the model-predicted trend for the fractional change in the relative 3.3 μm feature brightness (black line).
Download figure:
Standard image High-resolution image

5.1.2. Photodestruction Model

We simulate the effects of photodestruction of small PAHs as the average hardness of the radiation field increases inversely with metallicity by imposing a sliding minimum grain size below which all PAHs are destroyed. In this scenario, q_PAH and ΣPAH/TIR decrease because small PAHs are removed from the population starting from the smallest grains up. We model this with a minimum grain size cutoff, a_min, which rises as metallicity drops. The functional form of $\mathrm{log}({a}_{\min }(Z))$ is a power law, with index determined by minimizing χ² between q_PAH from the model GSDs and the observationally determined q_PAH values. The model is summarized in Figure 10. The best-fit correlation for a_min(Z) is

$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{{a}_{\min }}{4\mathring{\rm A} }\right)\,=\,\left\{\begin{array}{ll}0 & Z\gt {Z}_{\mathrm{th}}\\ 0.52\times \left[1-{\left(\displaystyle \frac{Z}{\,{Z}_{\mathrm{th}}}\right)}^{3.2}\right] & Z\leqslant {Z}_{\mathrm{th}}\end{array}\right.\end{eqnarray} \tag{ 5 }$

and is shown in the top-right panel of Figure 10. We find this model can be fit to the q_PAH trend and that it coarsely reproduces the ΣPAH/TIR decline. The resulting band ratio trends, however, are a poor match for the observations. Since q_PAH is constant at metallicities higher than Z_th, the modeled band ratios likewise do not change. We found this model does not explain many of the observed trends. The resulting shift in the PAH band ratios with metallicity is for many features the opposite of the observed trend, as can be expected for a scenario that preferentially weakens short wavelength bands. For instance, the bottom panel of Figure 10 shows that modeled PAH 6.2/ΣPAH decreases but the observed trend increases and PAH 17/ΣPAH increases dramatically in the model while the observed trend drops steeply. The model does reproduce the drop in the 7.7 μm feature below Z_th; however it does not explain the observed rise in this feature from Z_⊙ to Z_th. Even a combination of the radiation hardness and photodestruction models cannot produce the observed shift of power from long- to short-wavelength PAH features because the small PAHs that emit most effectively at shorter wavelengths are being preferentially destroyed faster than the radiation hardness is increasing.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Grain photodestruction model summary. Top left: the size distribution of PAH-emitting grains varies as a function of metallicity, with the lower size cutoff a_min increasing as metallicity drops below a specific threshold, at a rate (top right) that reproduces the q_PAH-Z relationship. The dashed blue line indicates grains with N_C = 10³ so at this a_min, q_PAH = 0. Bottom: comparisons between six observed PZR trends (black) and the photodestruction model (red). The dashed vertical line indicates the threshold metallicity below which a_min begins increasing. Seven metallicity bins indicate the median of each trend in cyan, with detections below 3σ set to zero.
Download figure:
Standard image High-resolution image

5.1.3. Inhibited Growth Model

We simulate the effects of inhibited PAH growth by reducing the abundance and fiducial size of PAHs. The fiducial size parameter a₀₁ characterizes the location of the first peak in the bimodal log-normal GSD of carbonaceous dust grains seen in Figure 8. To reproduce the observed q_PAH, we also use a multiplicative factor B_r that rescales the entire GSD as metallicity falls.

The functional form of a₀₁(Z) is assumed to be a power law: a₀₁ = 5 Å (Z/Z_⊙)^β, where the value of a₀₁ at Z_⊙ is chosen to match the "large" end of sizes considered in D21 and β > 0. We then explore the parameter space of β ≤ 6 in steps of 0.5 by calculating the resulting χ² between the the model and six primary observed PZR trends (i.e., q_PAH, ΣPAH/TIR, PAH 6.2/ΣPAH, PAH 7.7/ΣPAH, PAH 11.3/ΣPAH, and PAH 17/ΣPAH). We assume a similar power-law form for B_r(Z) but with a different exponent, constrained by minimizing χ² between the model and q_PAH trends.

Figure 11 summarizes the inhibited growth model that most closely matches observations. Since q_PAH is held constant as metallicity drops from Z_⊙ to Z_th, the amplitude of the GSD stays approximately constant until metallicity drops below Z_th. The best value for β is 1, with negligible improvements from increasing to 1.5 or 2; future work with more data will allow us to better constrain this parameter.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Inhibited grain growth model summary, as in Figure 10. Top left: GSDs as a function of metallicity. Top right: mapping between average PAH size a₀₁ (black) and GSD rescale factor B_r (pink) and metallicity. Bottom: comparisons between six observed PZR trends (black) and our best inhibited growth model (red).
Download figure:
Standard image High-resolution image

Since B_r is constrained based on the q_PAH–metallicity trend, the functional form of B_r is equivalent to a fit to the q_PAH trend. Above Z_th, B_r increases slightly with metallicity to keep q_PAH constant while a₀₁ shifts to larger sizes. The variable q_PAH has a small dependence on a₀₁ because it is defined as the integral of the GSD for grains with N_C < 10³, so as a₀₁ is increased, more of the distribution falls above this range and q_PAH decreases slightly. Below Z_th, q_PAH decreases as a power law with index 2.1 and B_r has this same relation normalized to one at Z_th.

With a₀₁ decreasing linearly as metallicity drops, the relative changes in each major PAH feature are reproduced to first order (see Figures 11 and 12). Specifically, PAH 17/ΣPAH drops by 50% at 0.4 Z_⊙ relative to Z_⊙, PAH 11.3/ΣPAH drops by ∼20%, PAH 7.7/ΣPAH drops by ∼5%, and PAH 6.2/ΣPAH rises by ∼50%. There are residual differences between the best model and the observations. For example, it is not possible to match the strong increase in PAH 6.2/ΣPAH of almost 75% between Z_⊙ and 0.4 Z_⊙ with any value of β. Similarly, we cannot reproduce the ∼30% drop in PAH 7.7/ΣPAH after a 10% rise from Z_⊙ to Z_th. The model does show this rise and fall in PAH 7.7/ΣPAH, but the amplitude is much smaller. The photodestruction model results in a better match for this specific trend but fails significantly for other bands.

Figure 12. Refer to the following caption and surrounding text. — **Figure 12.** Fractional PAH power as a function of metallicity for our best inhibited growth model (Section 5.1.3). The data bins from Figure 7 are shown. The gray lines indicate extrapolations of the model to lower metallicities than are probed by our observations. Also shown is the model-predicted trend for the fractional change in the relative 3.3 μm feature brightness.
Download figure:
Standard image High-resolution image

The inhibited growth model does an excellent job reproducing the observed trend between ΣPAH/TIR and metallicity. This trend is not explicitly fit; the modeled trend is a result of employing the factor B_r to best match the q_PAH trend and fitting a₀₁(Z) to best match the ratios PAH/ΣPAH. While q_PAH is constant as metallicity drops from Z_⊙ to Z_th, ΣPAH/TIR rises because the average PAH size a₀₁ decreases and smaller PAHs emit more per unit mass (D21). Below Z_th, q_PAH begins to drop as well, and this causes ΣPAH/TIR to fall faster than the shifting of a₀₁ causes it to rise.

6. Discussion

6.1. Inferences from Observational Trends

As metallicity decreases from Z_⊙ to 2/3 Z_⊙, ΣPAH/TIR remains fairly constant at about 15% with a scatter of about 2% (see Figure 4). As metallicity decreases below 2/3 Z_⊙, ΣPAH/TIR begins to drop following a power law of index 2.6. The slope of this power law is steeper than observed in previous works. Rémy-Ruyer et al. (2015) found a slope of about 1.3 between their f_PAH and oxygen abundance for galaxies in the Dwarf Galaxy Survey and part of the KINGFISH sample (Kennicutt et al. 2011; Madden et al. 2013). The difference may be a result of sample variations since Rémy-Ruyer et al. (2015) compares photometric integrated galaxy measurements for the PAH fraction, while in this work we study the PAH fraction from spectroscopy within individual galaxies. Future studies in this series with a significantly larger sample will allow us to place better constraints on the power-law index and scatter of the ΣPAH/TIR-metallicity relationship.

Some of the H ii regions are notably offset below the overall ΣPAH/TIR trend with metallicity. Since these regions have a lower ΣPAH/TIR than their metallicity would otherwise imply, this may be evidence of a local PAH luminosity deficit inside ionized gas within H ii regions. Previous works have proposed that the hard radiation environment near young O and B stars destroys small dust grains (Povich et al. 2007; Chastenet et al. 2019; Egorov et al. 2023). In the UV-illuminated neutral and molecular gas of the photodissociation region surrounding the H ii region, the small grains survive and have a steady source of photons for excitation. At the ∼300 pc scales probed by the apertures used in this work, our H ii region pointings contain a mix of both of these environments, resulting in a moderately lower ΣPAH/TIR in some cases. However, intrinsic features of the PAH spectrum (i.e., band-to-total and band-to-band ratios) do not appear to be distinct for H ii and non-H ii regions at a given metallicity, consistent with the mild effects of ISRF intensity and hardness (see Figures 4 and 9, respectively). This could indicate that the PAH population remains fairly uniform outside H ii regions but that there is a deficit of all PAHs relative to larger grains in the ionized gas itself (e.g., Sutter et al. 2024).

The modeled quantity q_PAH is typically derived from the brightness at IRAC 8 μm, which only traces the 7.7 μm PAH feature. In this work we have found that all PAH feature strengths decrease rapidly relative to TIR when metallicity falls below Z_th but that some features drop faster than others. Specifically, relative to the total dust luminosity, the 6.2 μm feature decreases more slowly than the 7.7 μm feature and the 17 μm feature decreases significantly faster as metallicity drops, implying that the situation is more complex than an overall decrease in PAH abundance (q_PAH). The assumption that the 7.7 μm feature is representative of all PAH emission underlies many past interpretations based on 8 μm photometry. The common explanation for the decrease in 8 μm photometry in such studies is photodestruction of the smallest PAHs. However, our observations show an increase of 6.2 μm/ΣPAH at low metallicity (see Figure 7). The 6.2 μm feature is more dependent on the presence of small PAHs than 7.7 μm, so the increase in 6.2 μm/ΣPAH is evidence against this photodestruction explanation for the PZR trends.

The shift in power between PAH features provides a good diagnostic of what is happening to the PAHs as metallicity changes. We observed a general trend of power shifting from longer-wavelength features to shorter-wavelength features. This occurs immediately as metallicity decreases below Z_⊙, despite the fact that q_PAH and ΣPAH/TIR remain approximately constant (with bulge ISRF correction applied; see Section 4.1.1). Previous works have observed similar trends. Smith et al. (2007) found the ratio PAH 11.3/17 μm increases from ∼1.2 to ∼2.2 as metallicity decreases from 0.8 Z_⊙ to 0.45 Z_⊙. We find a more significant increase in PAH 11.3/17 μm over this range, from ∼1.5 to ∼10.1 (see top panel of Figure 6). Note the average at 0.45 Z_⊙ is strongly affected by the points from NGC 2403, while the points from M101 show a lower average of ∼2.5 that is more in line with the Smith et al. (2007) result.

Another nearby low-metallicity galaxy with extensive Spitzer-IRS spectroscopy is the SMC (∼0.2 Z_⊙; Toribio San Cipriano et al. 2017). Measurements of the PAH fraction and band strengths in the SMC are in good agreement with the trends predicted by our inhibited growth model at that metallicity and do not match predictions with the photodestruction or ISRF hardness models. In particular, q_PAH for the SMC is ≲1% (Sandstrom et al. 2010; Chastenet et al. 2019). The band strengths in the SMC show a similar shift of PAH power between features, namely, higher PAH 6.2/ΣPAH and PAH 11.3/ΣPAH and lower PAH 7.7/ΣPAH compared to higher-metallicity regions (Sandstrom et al. 2012). The SMC regions also showed notably low PAH 17/ΣPAH compared to the SINGS average. These trends are all in line with a scenario where inhibited growth leads to smaller PAHs and lower PAH fractions in the low-metallicity SMC.

In addition to metallicity, many physical properties vary with radius in spiral galaxies. Of relevance to PAH emission, changes in the radiation environment including the intensity and average photon hardness are expected. As found in Figure 5 and Section 5.1.1, these properties are unlikely to drive the observed trends. Other potentially radially varying properties we do not model that could affect PAH emission include the phase balance of dense and neutral gas, gross differences in the optical properties of small carbon grains, varying carbon depletion fractions, and the metal content of the gas that exists in the dense phase.

6.2. Modeling Interpretation

We tested three physically motivated explanations for the observed PZR trends using D21 models: ISRF hardness, photodestruction, and inhibited grain growth. A summary of how each of these models compares to the observed PAH band ratios is shown in Figure 13. The ISRF hardness model quantified the effect of changing only the incident starlight spectrum with a fixed PAH population. The photodestruction model tested the most common explanation for the PZR: selective destruction of small grains. The inhibited growth model approximated the effects of reduced accretion of carbon onto grains by decreasing the average size and overall abundance of PAHs.

Figure 13. Refer to the following caption and surrounding text. — **Figure 13.** Summary of observed and modeled PAH/ΣPAH differentials. Top: percentage change in each major PAH feature relative to the sum of all PAH emission, as metallicity decreases from Z_⊙ to Z_th. Observations are in black, and the three models considered appear as colored curves. Bottom: same as above but for metallicity decreasing from Z_th to 0.45 Z_⊙. Note that the models differ most in their predictions for the 3.3 μm feature, not included in the current data set.
Download figure:
Standard image High-resolution image

The ISRF hardness model was intended to simulate a scenario where the PAH GSD remains unchanged and the PZR trends arise purely as a result of the increasing radiation hardness as metallicity decreases. This model resulted in the correct band ratio trends, but the magnitude of these trends is too low by a factor of 2–3. Furthermore, the adopted model explores the extremes of such a hardness trend, varying the incident starlight spectrum from that of a 3 Z_⊙ old bulge down to a Z_⊙/20 starburst. We conclude that starlight hardness alone cannot account for the varying shape of the PAH spectrum and that the small GSD must also be changing.

If, however, the ratio of far-UV-to-optical absorption cross sections of the PAH grains were substantially increased, then such a model could better match the observations as grain heating rate would exhibit increased sensitivity to the starlight hardness. Donnelly et al. (2024) also concluded that such a change in the relative absorption cross sections would result in better agreement between the D21 models and PAH emission observed near a low-luminosity active galactic nucleus. Regardless, the radiation hardness model is not capable of reproducing the q_PAH and ΣPAH/TIR trends as a function of metallicity. As the average exciting photon energy increases, it is not possible for either q_PAH or ΣPAH/TIR to decrease. In fact, without changes in the grain population, relative PAH power is expected to increase with hardness, as every PAH grain reaches higher stochastic temperatures.

The photodestruction scenario we modeled in Section 5.1.2 is commonly invoked to explain the deficit of PAH emission relative to total dust luminosity as a function of decreasing metallicity (Madden et al. 2006; Gordon et al. 2008). To reproduce the steep q_PAH decline, the photodestruction model requires increasing the grain cutoff size a_min up to ∼400 N_C to match the q_PAH drop from its Z_⊙ value of ∼3.9% to below 0.5% at 0.2 Z_⊙. Based on theoretical work, it is highly unlikely for a majority of PAHs smaller than 400 N_C to be photodissociated in ISM at low metallicity. Guhathakurta & Draine (1989) found that PAHs smaller than ∼ 23 N_C cannot survive in the ISM under typical solar-neighborhood radiation conditions. In Guhathakurta & Draine (1989) the significantly harder radiation field of a B-type star is considered and is found to increase the minimum grain size to a_min ∼ 35 N_C. This is an order of magnitude too low to explain the observed trend in q_PAH with photodestruction alone.

Because of this and the strong disagreement in band ratio behavior, we conclude that small grain photodestruction likely does not play a dominant role in the PZR. We do note that this is based on PAH emission trends at wavelengths longer than 5 μm. It remains possible that a modest amount of photodestruction may affect the very smallest PAHs traced by the 3.3 μm feature, which the photodestruction model predicts will drop precipitously toward low metallicity.

The model that best matches the observed PAH–metallicity trends simulates inhibited grain growth. We implemented this scenario by shifting the average PAH size smaller and reducing the total mass of all carbon grains as metallicity decreases. If individual carbon atoms are accreted and grow on the surfaces of larger grains (bottom-up formation), their reduced availability and longer accretion timescales at lower metallicity could explain the overall shift toward smaller grain sizes, although the increased intensity and hardness of radiation could also play a role in inhibiting growth.

While the metallicity dependence of the gas-grain accretion timescale supports this interpretation, any physical mechanism that substantially shifts the PAH GSD to smaller sizes would produce similar effects on the PAH spectrum. For instance, grain sizes injected by carbon-rich AGB stars could vary with initial stellar metallicity (Otsuka et al. 2014), although the nature of any such dependency is not well established. Another possible mechanism is selective processing of the largest PAH grains. However, photoprocessing is believed to have the greatest impact on smaller PAHs, and grain erosion from selective sputtering of all large grains would serve to increase the fractional PAH luminosity at long wavelengths, which does not match the observations.

Despite the successes of the inhibited growth model, several trends in the observations remain unaccounted for. These may be a result of simplifying assumptions in the D21 models. We have assumed that the radiation field intensity remains fixed at U = 1. This is the typical value of U extracted from the SED modeling maps with our matched apertures, and there is little effect on the shape or relative strength of the PAH spectrum at such low intensities. On smaller scales (e.g., approaching individual H ii regions), however this assumption may break down and larger values of U that could lead to multiphoton heating may be present.

Similarly, aside from the radiation hardness model, we have assumed the mMMP radiation field of the solar neighborhood applies to all regions in each galaxy. We have no direct measure of the shape of the incident SED that illuminates the dust, but the radiation hardness model (Section 5.1.1) shows that even shifting from an M31 bulge SED to a Z_⊙/20 10 Myr SED cannot reproduce the observed band ratio variations as a function of metallicity. Since our data do not span such a wide range in metallicity, the local radiation fields are unlikely to be changing this rapidly. However, previous works have found evidence that band ratio variations are correlated with radiation field hardness in nearby galaxies (Baron et al. 2024).

All three models assume unchanging vibrational properties of PAH grain material; if the chemical structure of PAH grains change as metallicity declines—for example, losing the ability to emit effectively at 17 μm—this could offer an alternative explanation for the success of the inhibited growth model.

Finally, we have also made the assumption that the ionization fraction as a function of grain size remains fixed at the standard value. This was done primarily to reduce the number of free parameters. However, the overall success of our inhibited growth model implies the standard ionization fraction as a function of size is reasonable. Since the ionization fraction is lower for smaller grains and the average grain size in our inhibited growth model decreases, the model naturally results in reduced effective ionization fraction. This is consistent with the observations of the PAH 7.7/11.3 μm ratio shown in Figure 6.

The strong agreement between the inhibited growth model and the observations emphasizes the role of accretion and regrowth in maintaining the small grain population in the ISM. At low metallicities, these processes become less efficient as the carbon abundance and average size and number of large grains drop (Zhukovska et al. 2008). This suggests that the majority of small grains are formed by accretion in the denser phases of the ISM and are not directly stellar in origin.

6.3. Future Work

The spectral coverage of Spitzer did not include the aromatic feature at 3.3 μm. Theoretical and laboratory results show this feature originates almost exclusively from the smallest PAH molecules (<100 N_C; Maragkoudakis et al. 2020; D21). If these small grains are significantly affected by photodestruction, we expect that the 3.3 μm emission will become very weak at sufficiently low metallicity. And yet, the low-metallicity galaxy IIZw40 (∼20% Z_⊙) has been found to have a relatively high fractional 3.3 μm luminosity (∼3%; Lai et al. 2020).

We have also seen little evidence in this work that small PAHs are preferentially destroyed at low metallicity. In fact, the success of our inhibited growth model suggests that even as the total PAH grain population drops in abundance, smaller PAHs rise in relative importance. However, the 3.3 μm feature is far more sensitive to the abundance of PAHs with N_C ≪ 100 than the 6.2 or 7.7 μm features (Hensley & Draine 2023), so our conclusions in this work cannot rule out photodestruction thresholds at smaller grain sizes.

If our inhibited grain growth model holds true, we predict the 3.3 μm fraction should carry a higher fraction of the total emission from PAH bands as metallicity decreases. The best inhibited growth model suggests the 3.3 μm feature contributes only 2% of all PAH band emission at solar metallicity, but extrapolating to 0.2 Z_⊙, we predict 3.3 μm could contribute about 10% of all remaining PAH band emission. As shown in Figure 13, the relative behavior of the 3.3 μm PAH feature has the strongest variation between models. Future studies on the PZR with 3.3 μm using JWST can disentangle which model or combination of models applies to the smallest PAHs with N_C ≪ 100.

Figure 4 shows ΣPAH/TIR increases from 10% to 16% as metallicity drops from Z_⊙ to Z_th. We applied a correction for this in Section 4.1.1 before performing the fit under the assumption that the ISRF is softer in the central region of M101. However, we also found our inhibited growth model naturally results in an excellent match to the observed ΣPAH/TIR trend as the average grain size shifts larger at higher metallicity. Future work with high-metallicity regions dominated by young stars can explore this further and determine whether the decrease in ΣPAH/TIR is a result of softer radiation environment or grain growth in high-metallicity ISM. If the trend is truly metallicity dependent, this may suggest that the threshold metallicity 2/3 Z_⊙ forms a modest local maximum in PAH brightness relative to all other dust.

The high luminosity of PAHs, dropping steeply and changing shape below ∼2/3 Z_⊙, may indicate that PAH emission could serve as a powerful probe of metal content in galaxies. JWST has been able to directly detect 3.3 μm PAH emission at z ∼ 4.1 (Spilker et al. 2023), and evidence from UV attenuation at z ∼ 7 may indicate the presence of large quantities of carbonaceous grains at very early times (Witstok et al. 2023). A future far-IR observatory with spectroscopic capabilities could offer an effective way to recover dust conditions and metal content in the early Universe (Moullet et al. 2023).

7. Conclusions

Spitzer-IRS spectroscopic radial strip maps across three nearby spiral galaxies with well-characterized metallicity gradients were used to study the deficit of PAH emission relative to total dust luminosity as metallicity declines. We found ΣPAH/TIR remains relatively constant between Z_⊙ and 2/3 Z_⊙. Below this threshold, ΣPAH/TIR drops smoothly and rapidly with metallicity. We studied how the PAH emission spectrum changes in strength and form to diagnose the physical origins of the PAH–metallicity relationship. Building on the latest models from D21, we reproduced these trends with a model of declining grain size and overall abundance designed to simulate the effects of inhibited PAH growth. The main findings of this work are as follows:

1.
PAH emission constitutes about 15% of all infrared emission in high-metallicity regions of normal star-forming galaxies. Below about 2/3 Z_⊙, the power emitted in the PAH features begins to drop smoothly with metallicity, approximately as a power law with slope ∼2.6.
2.
As its overall luminosity declines with lower metallicity, PAH feature power shifts toward shorter wavelengths, with the 17 μm band declining particularly rapidly. The dominant 7.7 μm band rises, but then it drops below the threshold metallicity, potentially due to the competing effects of smaller average grain size and the implicit changes in PAH ionization balance that naturally follow.
3.
Models based on the commonly assumed scenario of small grain photodestruction—when tuned to match measured trends in q_PAH—fail to reproduce observed changes in the PAH spectrum, predicting roughly the opposite band ratio trends as observed. We conclude that photodestruction cannot be the dominant effect driving the decline in PAH power at low metallicity.
4.
Although the harder radiation at low metallicity does shift PAH power to shorter wavelengths, the effects are modest, and modeling suggests the observed trends cannot arise solely from this increased radiation field hardness. Strong changes to the underlying carbonaceous GSD are needed to reproduce these trends.
5.
A model of inhibited grain growth in which the average size of PAH grains declines while the overall PAH abundance drops can explain most details of the varying PAH spectrum toward low metallicity. This may be evidence of the importance of dust regrowth as grains cycle through dense phases of the ISM.
6.
Our best inhibited growth model predicts the 3.3 μm PAH feature strength relative to the total PAH luminosity will increase from 2% to 10% as metallicity decreases from Z_⊙ to 0.2 Z_⊙, although moderate levels of photodestruction that do not substantially impact longer-wavelength bands could limit this value.

Acknowledgments

C.W. thanks R. Chandar's University of Toledo Paper Writing course (2022–2023) for constructive feedback. We thank K. Gordon, C. Engelbracht, E. Tarantino, D. Narayanan, E. Peeters, L. Hunt, and A. Li for helpful discussions over the years on the PAH–metallicity problem and D. Berg for insights on metallicity gradients. J.D.S. thanks the Max Plank Institut für Astronomie for hosting numerous scientific visits that contributed to the development of this work. We gratefully acknowledge support from NASA ADAP grant No. 80NSSC21K0851. This research has made use of NASA's Astrophysics Data System and Tiny-Tim/Spitzer, developed by John Krist for the Spitzer Science Center. This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA.