Estimating Molecular Gas Content in Galaxies from PAH Emission

Emission from polycyclic aromatic hydrocarbons (PAHs), a commonly used indicator of star formation activity in galaxies, also has the potential to serve as an effective empirical tracer of molecular gas. We use a sample of 19 nearby galaxies with spatially resolved mid-infrared Spitzer spectroscopy, multi-wavelength optical and mid-infrared imaging, and millimeter interferometric CO(1-0) maps to investigate the feasibility of using PAH emission as an empirical proxy to estimate molecular gas mass. PAH emission correlates strongly with CO emission on sub-kpc scales over the diverse environments probed by our sample of star-forming galaxies and low-luminosity active galactic nuclei. The tight observed correlation, likely a consequence of photoelectronic heating of the diffuse interstellar gas by the PAHs, permits us to derive an empirical calibration to estimate molecular gas mass from the luminosity of PAH emission that has a total scatter of only {\sim} 0.2 - 0.25 dex. Mid-infrared bands sensitive to PAH emission (e.g., the Spitzer/IRAC4 and WISE/W3 filters) can also be used as a highly effective substitute for this purpose.


INTRODUCTION
The formation of stars from dense molecular clouds in the interstellar medium constitutes an essential process in the lifecycle of galaxies. The physical connection between star formation and the interstellar medium can be explored through many observational frameworks, among them the tight correlation between SFR surface density (Σ SFR ) and cold gas surface density (Σ gas ). Originally proposed by Schmidt (1959) and subsequently extended and refined most notably by Kennicutt (1998), the Kennicutt-Schmidt law stipulates that star formation in galaxies conforms to a power-law relation of the form Σ SFR ∝ Σ N gas , where N ≈ 1.4 if both atomic (H I) and molecular hydrogen (H 2 ) are included, or N ≈ 1.0 if only H 2 is included (e.g., Bigiel et al. 2008;Leroy et al. 2013). To fully understand the physical processes that regulate galaxy formation and evolution, large galaxy samples with reliable measurements of SFR and gas content are indispensable.
While a broad range of SFR estimators has been extensively developed and applied to normal star-forming galaxies (see reviews in Kennicutt &Evans 2012 andMadau &Dickinson 2014) and galaxies hosting active galactic nuclei (AGNs; e.g., , the measurement of molecular gas content is quite limited owing to the difficulties of direct H 2 detection. In practice, most of what we know about the molecular gas in galaxies comes from observations of "tracer" species. The most commonly employed tracer of H 2 is CO(1-0) emission (e.g., Young & Scoville 1991;Solomon & Vanden Bout 2005; see review in Saintonge & Catinella 2022), but in recent years increasing attention has been devoted to H 2 estimators based on dust emission (e.g., Magdis et al. 2011;Scoville et al. 2014;Shangguan et al. 2018) and absorption (Concas & Popesso 2019;Yesuf & Ho 2019Barrera-Ballesteros et al. 2020). The CO-based method relies on an uncertain CO-to-H 2 conversion factor (Blanc et al. 2013;Bolatto et al. 2013), and a metallicity-dependent dust-to-gas mass ratio needs to be invoked to translate the dust continuum emission to molecular gas mass (Galametz et al. 2011;Rémy-Ruyer et al. 2014).
PAH emission, which arises primarily from the warm layer of photodissociation regions in proximity to H II regions, is closely associated with both dense, star-forming gas and more diffuse interstellar clouds. Within a photodissociation region, most of the incident ultraviolet radiation from the H II region is absorbed by PAH molecules and dust grains of size somewhat larger than the PAHs (Bakes & Tielens 1994;Wolfire & Kaufman 2011). After absorbing an ionizing ultraviolet photon, a PAH molecule mostly undergoes radiative deexcitation through IR emission, which accounts for up to 20% of the total IR power of galaxies (Smith et al. 2007;Diamond-Stanic & Rieke 2010;Xie et al. 2018). PAH molecules also can deexcite through photodissociation or ionization (Allamandola et al. 1989;Allain et al. 1996), in the process converting 0.1%-1% of the absorbed energy to energetic photoelectrons that dominate the heating of photodissociation regions as well as the diffuse interstellar clouds (Bakes & Tielens 1994;Weingartner & Draine 2001). These sources of heating in the photodissociation region are balanced by cooling through collisionally excited fine-structure lines of [C II] 158 µm and [O I] 63 µm, as well as collisionally excited rotational-vibrational lines of H 2 and rotational lines of CO and other abundant molecules (Tielens & Hollenbach 1985;Wolfire et al. 1989;Meijerink & Spaans 2005; see review by Hollenbach & Tielens 1997).
In addition to the above coupling between heating and cooling processes, Bendo et al. (2008Bendo et al. ( , 2010Bendo et al. ( , 2020 proposed, on the basis of the spatial coincidence between PAH and cold dust emission observed within galaxies, that PAH molecules can be mixed to varying degrees with the diffuse cold dust, and hence cold gas. Theoretical calculations and observations indicate that not only ionizing ultraviolet photons from newborn stars but also non-ionizing radiation from evolved stars contribute to the excitation of PAH molecules (e.g., Draine & Li 2001;Sellgren 2001;Haas et al. 2002;Li & Draine 2002;Zhang & Ho 2022). Under such circumstances, the interstellar radiation field simultaneously heats the PAHs and the diffuse interstellar dust, and, so too, the gas. The intimate mixture of PAHs and cold dust ensures a tight physical coupling between PAHs and the gas component in galaxies. Moreover, PAH emission, as an effective indicator of star formation activity in diverse environments (e.g., Calzetti et al. 2007;Farrah et al. 2007;Pope et al. 2008;Treyer et al. 2010;Shipley et al. 2016;Maragkoudakis et al. 2018;Xie & Ho 2019), may be linked especially closely to the cold molecular gas from which new stars form.
Within this backdrop, we wish to explore the degree to which PAH emission can be used as a proxy to estimate molecular gas content in galaxies. This question is especially timely in light of the enormous improvements in MIR spectroscopy afforded by the successful operation of the James Webb Space Telescope (Rigby et al. 2022). Indeed, previous studies using the Spitzer Space Telescope (Houck et al. 2004) already have suggested that in nearby galaxies PAH emission correlates with molecular gas (e.g., Regan et al. 2006;Pope et al. 2013;Cortzen et al. 2019;Gao et al. 2019;Chown et al. 2021), although the correlation is non-linear (e.g., Regan et al. 2006;Bendo et al. 2008) and the two components appear to be spatially offset (e.g., Casasola et al. 2008;Schinnerer et al. 2013). These previous investigations were based on broad-band photometry using MIR filters that contain strong PAH features (e.g., IRAC4, which contains PAH 7.7 and 8.6 µm; W3, which contains PAH 8.6, 11.3, and 12.7 µm), which cannot separate the emission features from the underlying continuum, or they relied on spatially integrated spectroscopy with the Spitzer Infrared Spectrograph (IRS; Werner et al. 2004). Our study improves upon previous efforts by taking advantage of spatially resolved analysis of the PAH and molecular gas using a sample of nearby galaxies with available mapping-mode IRS observations as well as interferometric maps of the cold gas (CO and H I). We explore the relationship between PAH and CO emission, taking into consideration their mutual correlation with SFR, with the ultimate aim of offering a new, improved empirical predictor of molecular gas content based on the observed strength of PAH emission and star formation activity.

Sample and Data Sets
The sample is selected from the Spitzer Infrared Nearby Galaxies Survey 1 (SINGS; Kennicutt et al. 2003), whose mapping-mode IRS observations provide us with the requisite spatially resolved PAH data on 75 nearby star-forming galaxies (SFGs) and galaxies hosting low-luminosity AGNs. SINGS provides a rich repository of multi-wavelength optical and IR images useful for our analysis (SINGS Team 2020). The mapping-mode observations of SINGS cover the central ∼30 ×50 region of each galaxy over the wavelength range ∼ 5 − 38 µm with a spectral resolution of λ/∆λ ≈ 64 − 128, depending on the observing module of IRS. The point-spread function (PSF) of IRS increases with wavelength, from full-width at half-maximum (FWHM) ∼ 1. 5 − 3. 5 at ∼ 5 − 14 µm to FWHM ≈ 4. 5 − 9. 0 at ∼ 14 − 38 µm.
To evaluate the relative effectiveness of estimating PAH strength through MIR broad-band imaging, we analyze images acquired by SINGS using the Spitzer IRAC4 band (λ eff = 7.87 µm; Fazio et al. 2004), as well as images taken with the W3 band (λ eff = 11.56 µm) as part of the all-sky survey of the Wide-field Infrared Survey Explorer 2 (WISE; Wright et al. 2010; NASA/IPAC Infrared Science Archive 2020). We use continuum-subtracted narrow-band Hα images, corrected for contamination by [N II] λλ6548, 6584 in the Hα bandpass (Kennicutt et al. 2008), to trace the distribution of SFR. The IRAC4, W3, and Hα images have PSFs with FWHM ≈ 2 , 6. 5, and 2 , respectively. The galaxy M 51 deserves special mention. Although formally contained in SINGS, M 51 has IRS mapping-mode observations with more extensive spatial coverage (∼40 ×240 ) from a separate, dedicated observations (program 20138; PI: K. Sheth; IRSA 2022) 3 . These data were analyzed by Zhang et al. (2021) and will be used here.  To secure spatially resolved CO maps of comparable resolution, we cross-matched SINGS with the BIMA Survey Of Nearby Galaxies 4 (BIMA SONG; Helfer et al. 2003; NASA/IPAC Extragalactic Database 2019), which was carried out with the 10-element Berkeley-Illinois-Maryland Association (BIMA) millimeter interferometer (Welch et al. 1996). BIMA SONG systematically imaged the 3 mm CO(1-0) emission of the central ∼190 of 44 spiral galaxies at an hydrogen. Of the 34 nearby galaxies covered by THINGS, 13 are in common with our sample. Zhang et al. (2021Zhang et al. ( , 2022 devised an optimal strategy to extract from low-resolution IRS mapping-mode observation spatially resolved measurements of the integrated (5 − 20 µm) PAH emission, as well as the strengths of prominent individual PAH features. Where the IRS observations suffered from incomplete wavelength coverage, MIR photometry provided complementary coverage that still enabled effective spectral decomposition. Here we only give a very brief summary of the technical details.  After fitting and subtracting the residual background with a two-dimensional polynomial function, all the ancillary MIR and Hα images were convolved to a common angular resolution, which is constrained by the coarsest PSF of the IRS spectral data cube (FWHM ≈ 9. 0). We matched the resolution of the CO and H I maps in a similar manner, using the Gaussian2D function in astropy.modeling.models (Astropy Collaboration et al. 2013) to mimic their original elliptical synthesized beam and the create matching kernel function in photutils 6 to construct the convolution kernels. The convolution was performed using the convolve fft function in astropy.convolution. After convolving all the data to the same angular resolution, we reprojected them into the same coordinate frame as the PAH map to facilitate the subsequent pixel-based, spatially resolved analysis. The final 6 https://media.readthedocs.org/pdf/photutils/stable/photutils.pdf pixel size of the reprojected data was set to 10 to ensure that all the pixels within each galaxy are spatially independent; this size corresponds to a physical scale of ∼ 0.2 − 0.7 kpc for the distances of our sample (Table 1). Some pixels might require additional optimal binning to achieve the minimal signal-to-noise ratio required for robust PAH measurements (Zhang et al. 2021). The final number of spatially resolved spaxels for each galaxy ranged from 1 to 15, with the exception of M 51, which is ∼ 100.

Data Processing and PAH Measurement
The PAH strength was measured using the template-fitting technique of Xie et al. (2018) from complete (SL+LL) lowresolution IRS spectra. After masking out the isolated ionic emission lines, the IRS spectrum is fit with a multicomponent model consisting of a theoretical PAH template and three continuum components represented by modified blackbodies of different temperatures, all subject to attenuation by foreground dust extinction. The fitting is carried out with the Bayesian Markov chain Monte Carlo procedure emcee in the Python package, as described in Section 4.2 of Shang-guan et al. (2018). The median and standard deviation of the posterior distribution of each best-fit parameter are taken as the final estimate and corresponding uncertainty. Based on the best-fit model for each spaxel, we define the total PAH luminosity, L PAH , as the integral of the best-fit PAH template over the 5−20 µm region, following the convention of Zhang et al. (2021). Xie et al. (2018) showed that the measurements of individual PAH features based on the method used here are well consistent with measurements from other methods. The template-fitting method obtains ∼ 30% more total PAH flux compared to the PAHFIT (Smith et al. 2007) and CAFE (Marshall et al. 2007) methods because it better recovers the long-wavelength continuum associated with PAH skeletal vibrations.

The Correlation between Molecular Gas and PAHs
Previous work on PAH emission as a proxy for molecular gas was based on either global PAH measurements of galaxies (e.g., Pope et al. 2013;Cortzen et al. 2019;Gao et al. 2019) or spatially resolved measurements that use MIR broad-band photometry, which is sensitive to but does not give precise measurement of PAH emission (e.g., Regan et al. 2006;Bendo et al. 2010;Chown et al. 2021). Here we investigate, for the first time, whether PAH emission can trace molecular gas using accurate PAH measurements derived from spectral decomposition of spatially resolved MIR spectra (Section 2). We then used these improved mesurements to attempt to delineate the underlying mechanism for the tight correlation between PAH emission and molecular gas.
Following Solomon et al. (1997), we calculate the CO(1-0) luminosity, in units of erg s −1 , as where D L is the luminosity distance in cm, ν obs = ν rest /(1+ z) is the observed CO(1-0) line frequency, where ν rest = 1.153 × 10 11 Hz and z is the redshift, c is the speed of light in km s −1 , and S CO ∆V is the velocity-integrated flux in Jy km s −1 , which is converted from the reprojected moment 0 map in Jy km s −1 beam −1 by multiplying by 4ln2θ 2 πFWHM 2 , with θ the pixel size of the projected map and FWHM the beam size of the original map. The associated uncertainty is with the root-mean-square (rms) noise of the channel map in Jy, ∆v the channel width in km s −1 , and N the number of channels over which the line is integrated (Helfer et al. 2003). Note that only the spaxels with ∆ log L CO < 1 dex are used in this paper. The surface brightness of the CO emission is defined as Σ CO = L CO cos i/A, where i is the galaxy inclination angle and A is the spaxel area in pc 2 . Similar definitions apply to H I. We correct the Hα emission for Galactic extinction using the extinction curve of Cardelli et al. (1989) Schlafly & Finkbeiner (2011). We account for internal extinction following the prescription of Calzetti et al. (2007) for resolved H II regions, L Hα = L obs Hα + 0.031 L 24 µm , where L obs Hα is the Galactic extinctioncorrected Hα luminosity and L 24 µm is derived from MIPS24 images. Then, we adopt the prescription SFR = 5.3 × 10 −42 L Hα (Kennicutt & Evans 2012), assuming implicitly that it applies to galaxies on small ( 1 kpc) scales. As discussed by Kennicutt & Evans (2012), statistical approximations used to estimate SFRs may break down on small scales, but it is beyond the scope of this work to address this issue.
with a total scatter of t = 0.239 dex and an intrinsic scatter of i = 0.081 dex. The nearly linear correlation between CO and integrated PAH emission confirms the feasibility of using PAHs as a proxy for estimating molecular gas, even on sub-kpc scales. For completeness, Appendix A presents the results of substituting the spectroscopically derived integrated PAH emission with photometric estimates based on IRAC4 and W3 images. By sharp contrast, the H I line shows no relation with PAH emission whatsoever (Figure 2b). The lack of a connection between atomic hydrogen gas and PAHs may not be unexpected, given that the radial profile of H I differs from that of molecular hydrogen (e.g., Bigiel & Blitz 2012;Wang et al. 2014;Casasola et al. 2017). Spatially resolved studies of the Large Magellanic Cloud also reveal that CO correlates poorly with H I (Wong et al. 2009). What is the underlying physical mechanism behind the tight correlation between CO and PAH emission? Could it stem simply from their mutual connection with SFR? After all, SFR empirically correlates with PAH emission, and, at the same time, SFR is physically tied to the molecular gas via the Kennicutt-Schmidt law. Or, perhaps, the CO-PAH correlation reflects a common excitation mechanism (Section 1). Figure 3 shows the relationship between the surface brightness of the extinction-corrected Hα emission and that of the CO emission. As anticipated, the data confirm that the spatially resolved spaxels of the SINGS galaxies obey the Kennicutt-Schmidt relation. A linear regression gives with t = 0.362 dex and i = 0.220 dex. The slope of Equation 4 is consistent with a power-law index of N ≈ 1 for the Kennicutt-Schmidt relation involving only molecular gas, as obtained in previous studies of nearby galaxies spatially resolved on comparable scales (e.g., Bigiel et al. 2008;Leroy et al. 2013). At the same time, as is well-known, Hα emission, used here to indicate SFR, also correlates strongly with PAH emission (Figure 4; t = 0.208 dex, i = 0.175 dex): log Σ Hα = (1.030 ± 0.036)(log Σ PAH − 34.0) + (32.454 ± 0.053).
The closely similar slopes of Equations 3-5 do not offer much guidance as to whether the CO-PAH correlation is physical or merely an artifact of the SFR-CO and SFR-PAH relations. Nevertheless, two arguments can be advanced to suggest that molecular gas does have a direct link with PAHs. On the one hand, among the three empirical relations presented above, that involving CO and PAH (Equation 3) by far has the lowest formal intrinsic scatter (0.08 dex). We should be cautious, however, not to overinterpret this result, because the total scatter of all three correlations is likely dominated by many sources of systematic uncertainty, which are difficult to quantify. The exact scatter also depends on the statistical method used for the linear regression. On the other hand, the Spearman correlation between Σ CO and Σ PAH is highly significant, for all galaxies combined or for SFGs and AGNs separately (correlation coefficient ρ ≈ 0.80 − 0.90; p-value < 10 −20 ). Moreover, it still remains highly significant after taking into consideration the mutual correlation of these two variables with Σ Hα ; the partial rank correlation between Σ CO and Σ PAH while holding Σ Hα fixed drops only slightly (ρ ≈ 0.60; p-value < 10 −15 ). We suggest that the tight correlation between PAH and CO emission is physical and not an artifact of their mutual correlation with SFR. As Figure 3. The correlation between the surface brightness of extinction-corrected Hα emission and CO(1-0) emission. Red squares and blue circles correspond to spaxels from the central regions of AGNs and SFGs, respectively, while the corresponding symbols with black borders denote the nuclear (AGN-like) and offnuclear (SFG-like) spaxels from M 51. The black-dashed line indicates the best-fit linear regression (Eq. 4). mentioned in the Introduction, the photoelectronic effect of PAH molecules contributes to the heating of the gas (Wolfire et al. 1989;Stacey et al. 1991), and this might be the physical basis for the CO-PAH correlation. And while PAHs traditionally are viewed as concentrated toward photodissociation regions, they, in actuality, can be excited by the radiation field of the total stellar population and have a more extended spatial distribution, further enhancing their connection to the broader cold gas component (e.g., Jones et al. 2015;Bendo et al. 2020).
Compared to Equation 3, Equation 6 gives only a slightly reduced scatter ( t = 0.226 dex, i = 0.077 dex). Having argued that CO(1-0) emission is strongly correlated with PAH emission, and that their correlation may have a physical basis, we propose an empirical prescription to estimate molecular gas mass from the observed strength of PAH emission. We offer two prescriptions, one that relies solely on the strength of the PAH emission, and another that further introduces the SFR as a secondary variable. In this work, we use the extinction-corrected Hα emission to estimate SFR, but other estimators can be substituted. Adopting a canonical CO-to-H 2 conversion factor of α CO = 3.2 M pc −2 (K km s −1 ) −1 (this value should be multiplied by a factor of 1.36 to include helium; Bolatto et al. 2013), and converting Σ CO (erg s −1 pc −2 ) to velocity-integrated antenna temperature T int (K km s −1 ) using T int = 10 3 c 3 2kν 3 rest ΣCO 4π pc 2 , where k is Boltzmann's constant and pc is one parsec unit, we obtain, analogous to Equation 3, log M H2 = (1.029 ± 0.040)(log L PAH − 40.0) + (6.399 ± 0.040),

Estimating Molecular Gas Content from PAH Emission
where for which t = 0.226 dex and i = 0.080 dex ( Figure 5). Appendix A demonstrates that in the absence of spectroscopy even broad-band photometry using the Spitzer/IRAC4 and WISE/W3 filters can provide quite accurate predictors of CO emission. Table 3 summarizes the best-fit parameters for the M H2 calibration based on different measures of PAH strength, with and without including the effect of SFR as a second parameter.

SUMMARY
We study the correlation between PAH and CO(1-0) emission in galaxies, with the aim of elucidating the underlying physical connection between these two components of the interstellar medium and to explore the possibility of using PAH emission as a proxy for estimating molecular gas content. We use the method of Zhang et al. (2021) to extract robust, spatially resolved measurements of the total (5−20 µm) PAH emission by combining mapping-mode Spitzer/IRS observations with mid-IR images. We focus on the subset of 19 nearby galaxies that overlap between the SINGS and BIMA SONG projects, the former providing the Spitzer MIR spectroscopy and ancillary optical-MIR images necessary to measure PAH fluxes and SFRs, and the latter interferometric CO(1-0) maps of comparable resolution to derive molecular gas masses. H I maps of similar resolution from the THINGS collaboration are available for 13 of the 19 galaxies. The sample contains a mixture of star-forming galaxies and low-luminosity AGNs. The observations map the central ∼30 ×50 region of each galaxy, except for M 51, which has a larger coverage of ∼40 ×240 , at a spatial resolution of 10 , which corresponds to a physical scale of ∼ 0.2−0.7 kpc.
We confirm that PAH emission is tightly correlated with CO emission on sub-kpc scales, for both star-forming and active galaxies, as well as for different environments across a single galaxy. The same does not hold for H I gas. We argue that the correlation between CO and PAH emission is not an artifact of their mutual relation to SFR. Instead, we suggest that the CO-PAH correlation reflects an underlying physical connection between PAHs and the diffuse interstellar gas, likely tied to a common heating mechanism. We present an empirical calibration to estimate molecular gas mass from the luminosity of PAH emission that has a total scatter of only ∼ 0.2 − 0.25 dex. Including the SFR as a second parameter in principle can improve the accuracy of the molecular gas mass prediction, but the current sample is too limited to yield a firm conclusion on the merits of this approach. Calibrations using photometric estimates of PAH strength based on the Spitzer/IRAC4 and WISE/W3 filters give surprisingly similar results.  As has been investigated extensively, some broad-band MIR filters (e.g., Spitzer/IRAC4: Regan et al. 2006, Bendo et al. 2010; WISE/W3: Gao et al. 2019, Chown et al. 2021) are sensitive to the most prominent PAH features, offering a highly efficient alternative method to estimating PAH emission. These methods are attractive, in light of their applicability to large data sets available from large-area sky surveys. Figure A1 confirms that within our subsample of spatially resolved SINGS galaxies, the surface brightness of CO emission correlates strongly with the surface brightness measured in the IRAC4 and W3 bands. The best-fit regressions (Equations A1 and A2) have very similar slope and scatter compared to those obtained using spectroscopically derived PAH measurements (Figure 2a; Equation 3). Here, we define L band = νL ν (band), and we follow the prescription of Helou et al. (2004) to correct the IRAC4 band for stellar emission based on the IRAC1 map, with a scaling factor of 0.232. As a comparison, Figure A1b also presents the pixel-based calibration between the surface brightness of CO emission and that in the W3 band, as presented by Chown et al. (2021). Note that the calibration of Chown et al. (2021) is based on measurements of relatively lower surface brightness, as most of their resolved pixels are from the outer regions of galaxies.