The Calibration of Polycyclic Aromatic Hydrocarbon Dust Emission as a Star Formation Rate Indicator in the AKARI NEP Survey

From 2008-2016, we observed 19 slitmasks with Keck II/DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al., 2003) targeting the AKARI NEP-Deep field. When selecting targets for each slit in the later years, we prioritized mid-infrared AKARI/IRC sources with detected flux measurements at 9 and 18 $\mu$m, and assigned secondary priorities to high-redshift galaxies (photometric redshift $z_{phot}\sim 1$; Oi et al., 2014) and optical counterparts of Chandra X-ray sources (Krumpe et al., 2015). The spectrograph was configured with the 600ZD grating and GG495 order blocking filter (central wavelength $\lambda_{c}=7500$ Å) for the two slitmasks observed on August 24, 2014, and the 900ZD grating and GG455 order blocking filter ($\lambda_{c}=6800$ Å) for the remaining four slitmasks observed in 2015-2016. A typical slitmask contained $\sim 80$ slits that were 1^′′ in width. The spectral range covered approximately 4850 Å $\lesssim\lambda\lesssim 10050$ Å, depending on the location of the slit. The spectral resolution was $R\sim 2000$, which corresponds to full-width at half-maximum (FWHM) $\Delta\lambda\sim 3.75$ Å at the central wavelength $\lambda_{c}=7500$ Å, as confirmed by the observed widths of night-sky lines. Prior to observing each night, we acquired one set of Ne, Ar, Kr, and Xe arc lamp spectra and three flat field exposures of the detectors illuminated by an internal incandescent lamp. In addition, we observed standard star Feige 15 (RA $01^{h}49^{m}09.5^{s}$, Dec $+13^{\circ}33^{{}^{\prime}}11.4^{{}^{\prime\prime}}$) in long-slit mode for flux calibration. We reduced the raw data from 2014-2016 using the UCB DEEP2 spec2d reduction pipeline (Cooper et al., 2012; Newman et al., 2013), which processed the flatfields and calibration arcs, fit a polynomial wavelength solution to the pixel locations of arc lines, subtracted night-sky emission, and extracted one- and two-dimensional spectra. The extended flux along the slit was included, but most spectra were unresolved in the spatial direction.

The AKARI NEP collaboration observed eight DEIMOS slitmasks in July 2011, with the focus of targeting galaxy cluster candidates with AKARI/IRC MIR detections in the NEP-Deep field (PI: Tomotsugu Goto; Shogaki et al., 2018). 553 targets were selected based on Subaru/Suprime-Cam (Imai et al., 2007; Wada et al., 2008) R-band and Z-band brightness cuts such that $Z<23$ mag and $R>21.5$ mag. Slitmasks contained $\sim 50-80$ slits of 1^′′ width. The instrument configuration consisted of the 600ZD grating and GG455 order blocking filter, yielding a spectral resolution of $R\sim 2000$. A set of calibration images consisting of a Ne, Ar, Kr, Xe arc lamp image and five flat field lamp images was acquired for each slitmask. Bright spectrophotometric standard stars BD+$33^{\circ}2642$ and Feige 110 were observed for 90 seconds each in long-slit mode. Weather conditions were clear on both nights, with 0.6^′′ seeing. The data were reduced using the spec2d pipeline described above (Shogaki et al., 2018).

We then flux-calibrated the one-dimensional object spectra from 2011-2016 using standard IRAF routines as follows. The raw frames were overscan-trimmed and bias-subtracted using the mscred.ccdproc routine. A normalized flat field was created using mscred.flatcombine, noao.twodspec.longslit.response, and textttmscred.mscjoin, and the standard star spectra were flat-fielded with mscred.ccdproc. A wavelength-calibrated standard star image was produced using noao.twodspec.longslit.identify, reidentify, fitcoords, and transform on the arc images to determine the wavelength solution and noao.imred.kpnoslit.doslit on the standard star. Then, the sensitivity function was derived using noao.twodspec.longslit.standard, sensfunc, and calibrate on the standard star. Lastly, the sensitivity function was input into the routine noao.onedspec.calibrate to apply extinction corrections and flux calibrations to the reduced one-dimensional object spectra.

To estimate the photometric accuracy of the flux-calibrated spectra, we calculated the average flux density of the spectra within the CFHT/MegaCam r-band filter (5690-6930 Å):

where $F$ is the galaxy spectrum in erg/s/cm²/Å, $T$ is the filter response function, $\lambda$ is the wavelength in Å, and $c$ is the speed of light in Å/s. We compared our measurements to observed photometric CFHT r-band flux densities and applied a typical offset correction $\sim 0.1-0.25$ dex to account for slit losses. After applying the correction, the emission line fluxes have an uncertainty of $\sim 10$%.

The sample contains 68 H$\alpha$ and 259 [O II] emission line candidates with at least four MIR detections/upper limits in the 7-24 $\mu$m filters. The ranges in redshift, stellar mass, and TIR luminosity of the combined unique 296 sources are $0.05<z<1.36$, $2.1\times 10^{8}$ $M_{\odot}<M_{*}<9.3\times 10^{11}$ $M_{\odot}$, and $1.1\times 10^{9}$ $L_{\odot}<L(TIR)<2.2\times 10^{12}$ $L_{\odot}$. Table 1 summarizes the newly reduced and calibrated multi-slit spectroscopy from this work.

In 2017, we targeted higher redshift galaxies in the AKARI NEP-Deep field, with the Multi-Object Spectrometer For Infra-Red Exploration (MOSFIRE; McLean et al., 2010, 2012) on the Keck I telescope. Based on photometric redshifts provided by the AKARI-NEP team (Oi et al., 2014), we observed four slitmasks in the J-band (1.153–1.352 $\mu$m) and one slitmask in the Y-band (0.972–1.125 $\mu$m) to target H$\alpha$/[NII] in AKARI/IRC sources at $0.7<z<1.1$ and $0.5<z<0.8$, respectively. Each mask had $\sim 25$ slits with slit width of 0.7^′′, and resolution of $R\sim 3300$. Weather conditions were clear and dry on August 2nd and mostly clear with some high clouds on August 3rd. Calibration images consisted of Ne, Ar arc lamps (1.5 sec each) and internal flat fields (11 sec each). Slitmasks were observed in Multiple Correlated Double Sampling (MCDS) mode with 16 reads, and a dither position offset of $\pm$1.5^′′ was used. Standard stars FS 139 (RA $16^{h}33^{m}52.9^{s}$, Dec $+54^{\circ}28{{}^{\prime}}23{{}^{\prime\prime}}$) and HD 162208 were observed in long-slit mode.

We used the standard MOSFIRE Data Reduction Pipeline (DRP) to reduce the spectroscopic data. The pipeline creates a pixel flatfield, performs slit edge tracing, calibrates the wavelength scale, performs background subtraction, rectifies the slits, and extracts the 1D spectra. We flux-calibrated the 1D spectra by using the standard star. First we smoothed the standard star spectrum and interpolated over telluric absorption lines. Then we modeled the star as a blackbody based on its effective temperature, and normalized the curve to the star’s Ks-band magnitude from the Two Micron All Sky Survey (2MASS), factoring in the filter response function. The factor to convert from counts/second to erg/s/cm²/Å is then proportional to the ratio of the normalized blackbody function to the calibration star spectrum. We accounted for slit losses by comparing each galaxy’s average flux density to CFHT/WIRCAM J-band photometry.

Our MOSFIRE sample contains 14 H$\alpha$ and 9 [O II] emission-line candidates with at least four MIR detections/upper limits in the 7–24 $\mu$m filters. The ranges in redshift, stellar mass, and TIR luminosity of the combined unique 15 sources are $0.54<z<1.04$, $8.5\times 10^{9}$ $M_{\odot}<M_{*}<1.1\times 10^{11}$ $M_{\odot}$, and $5.5\times 10^{10}$ $L_{\odot}<L(TIR)<8.6\times 10^{11}$ $L_{\odot}$.

In addition to our 2011-2016 Keck II/DEIMOS observations, we obtained reduced DEIMOS spectra and emission line flux measurements from observations in semester 2008A through the AKARI-NEP collaboration (PI: Toshinobu Takagi). Targets included 242 AKARI-NEP mid-IR-selected sources with Subaru/Suprime-Cam R-band magnitude $\lesssim$ 24 mag. Five slitmasks were observed using the 600ZD grating and GG495 order blocking filter. Standard stars HZ 44 and Feige 110 were observed with a long slit (1.5^′′ slit width) using the same grating and filter setup as the targets. Data were reduced using the spec2d pipeline, and then flux-calibrated with standard IRAF routines. [O II]$\lambda\lambda 3726,3729$, H$\beta$, [O III]$\lambda 5007$, H$\alpha$, and [N II]$\lambda 6584$ emission lines were fit using Gaussian models with the IDL/MPFIT routine. Double Gaussians were used to fit the resolved [O II] doublet, H$\beta$ emission with underlying absorption, and H$\alpha$ with broad and narrow line components in AGN. The continuum was modeled with a linear fit. H$\alpha$ and [N II]$\lambda\lambda 6548,6584$ were fit simultaneously with three Gaussians. The Takagi et al. sample contains 48 H$\alpha$ and 126 [O II] emission line candidates with at least four MIR detections in the 7–24 $\mu$m filters. The ranges in redshift, stellar mass, and TIR luminosity of the combined unique 152 sources are $0.09<z<1.35$, $9.5\times 10^{8}$ $M_{\odot}<M_{*}<3.1\times 10^{11}$ $M_{\odot}$, and $8.3\times 10^{8}$ $L_{\odot}<L(TIR)<3.1\times 10^{12}$ $L_{\odot}$.

We also include spectroscopic redshifts and emission line fluxes from the Shim et al. (2013) sample, which consists of spectroscopy of hundreds of low-redshift galaxies and AGN ($z\lesssim 0.4$) in the AKARI NEP-Wide field. This provides rest-frame UV/optical emission line detections from MMT/Hectospec and WIYN/HYDRA multi-object spectrographs. The authors primarily selected targets with MIR detections at 11 $\mu$m and 15 $\mu$m. Secondary target selection criteria are summarized in their Table 1, and include AGN candidates, Herschel FIR sources, and PAH-luminous galaxies. The Shim et al. sample contains 498 H$\alpha$ and 571 [O II] emission line galaxieses with at least four MIR detections/upper limits in the 7–24 $\mu$m filters and secure redshifts (i.e., quality flag of 4). The ranges in redshift, stellar mass, and TIR luminosity of the combined unique 711 sources are $0.03<z<1.29$, $1.9\times 10^{8}$ $M_{\odot}<M_{*}<9.9\times 10^{11}$ $M_{\odot}$, and $4.3\times 10^{8}$ $L_{\odot}<L(TIR)<3.2\times 10^{12}$ $L_{\odot}$.

We supplement our high-redshift MOSFIRE data with near-infrared H$\alpha$ detections from Oi et al. (2017), who studied the mass-metallicity relation in luminous infrared galaxies at $z\sim 0.9$. Their sample consists of Subaru/FMOS H$\alpha$/[N II] multi-object spectroscopy in the J-long band (1.11–1.35 $\mu$m), of AKARI mid-infrared sources in the NEP-Deep field. Target galaxies were selected to have detections at 11 $\mu$m, 15 $\mu$m, and/or 18 $\mu$m, and estimated photometric redshifts $z_{phot}\sim 1$. H$\alpha$/[N II] fluxes are given for 28 secure H$\alpha$ emitters and H$\alpha$ fluxes are given for 36 “non-secure” H$\alpha$ detections (i.e., objects where a single emission line was observed and assumed to be H$\alpha$.) Some of the non-secure H$\alpha$ objects have since been re-observed in optical wavelengths by AKARI NEP collaborators using MMT/Hectospec, WIYN/Hydra, Keck II/DEIMOS, and GTC/OSIRIS, allowing us to confirm that 10/36 are indeed H$\alpha$ detections based on the detection of other strong emission lines at shorter wavelengths (e.g., [OII], H$\gamma$, H$\beta$, [OIII]). In addition, we find that the emission line from object 61016583, assumed by the Oi et al. to be H$\alpha$, is [O III]$\lambda$5007 at $z=1.357$. For reference, we list the updated secure H$\alpha$ objects in Table 9 in the Appendix. Our SFR calibration sample includes 37 FMOS objects with H$\alpha$ detections and at least four MIR detections/upper limits in the 7–24 $\mu$m filters. The ranges in redshift, stellar mass, and TIR luminosity are $0.70<z<1.03$, $2.1\times 10^{9}$ $M_{\odot}<M_{*}<1.5\times 10^{11}$ $M_{\odot}$, and $8.1\times 10^{10}$ $L_{\odot}<L(TIR)<1.5\times 10^{12}$ $L_{\odot}$.

The AKARI/IRC slitless SpectroscoPIC surveY (SPICY) obtained infrared spectra of all sources with 9 $\mu$m flux density brighter than 0.3 mJy in 14 NEP fields of $10{{}^{\prime}}\times 10{{}^{\prime}}$. The primary goal was to study PAH emission in galaxies at $z\lesssim 0.5$ (Ohyama et al., 2018) within the NEP-Deep and Wide fields. Low-resolution spectra ($R\sim 50$) were acquired with the short-wavelength MIR camera on the IRC, which covers a wavelength range of $5-13\ \mu$m. Ohyama et al. identified 48 galaxies with detectable PAH 6.2, 7.7, and 8.6 $\mu$m features. All 48 galaxies have AKARI/IRC photometry, with 83% having detections in the N2, N3, N4, S7, S9W, S11, L15, and L18W filters ($\sim 2-18\ \mu$m). Of the 48 PAH galaxies in the sample, 41 galaxies have optical spectroscopic detections and 33 galaxies have detected H$\alpha$ emission from Keck II/DEIMOS (this work) or MMT/Hectospec (PI: Ho Seong Hwang). Figure 2 illustrates typical SPICY spectra of star-forming galaxies from Ohyama et al. (2018). We use these 41 galaxies as our validation sample, to test the consistency and accuracy of our photometrically-derived PAH luminosity measurements. The ranges in redshift, stellar mass, and TIR luminosity of the 41 sources are $0.06<z<0.49$, $2\times 10^{9}$ $M_{\odot}<M_{*}<2\times 10^{11}$ $M_{\odot}$, and $2.1\times 10^{9}$ $L_{\odot}<L(TIR)<3.9\times 10^{11}$ $L_{\odot}$.

Figure 3 shows examples of the quality of our DEIMOS and MOSFIRE spectra. The top panel shows a MOSFIRE galaxy detected in the J-band that has red- and blue-shifted lines in H$\alpha$ and [N II] due to rotational kinematics. The second panel is an example of a typical star-forming galaxy observed with DEIMOS with H$\alpha$, [N II], and [S II] emission lines. The third panel shows a Seyfert galaxy with double-peaked emission in H$\beta$ and [O III]. The signal-to-noise is also high enough to resolve the doublet in H$\gamma$ (not shown). The fourth panel shows broad blue wings in [O III] emission, suggestive of outflows from an AGN.

We used the Specpro software (Masters & Capak, 2011) to visually inspect each spectrum and to determine an initial redshift guess based on a template of emission and absorption lines. This initial redshift served as the seed for our python routine that used non-linear least-squares fitting. We used Gaussian fits to measure redshifts, emission line fluxes, and equivalent widths for [O II]$\lambda\lambda$3726,3729, H$\gamma$, H$\beta$, [O III]$\lambda\lambda$4959,5007, H$\alpha$, [N II]$\lambda\lambda$6548,6584, and [S II]$\lambda\lambda$6716,6731. In general, a linear local continuum was added. For H$\gamma$ and H$\beta$, the continuum was modeled with a Voigt profile in order to account for underlying stellar absorption.

For multi-component Gaussian fits, the following constraints were used in order to limit the number of free parameters. When resolved, the $[$O II$]$$\lambda\lambda$3726,3729 doublet was fit with a double Gaussian, with the central wavelength and width of [O II]$\lambda$3729 fixed relative to [O II]$\lambda$3726. The $[$O III$]\lambda\lambda$4959,5007 lines were fit with a double Gaussian such that the width of [O III]$\lambda$4959 was fixed to that of the stronger line $\lambda$5007, the central wavelength was fixed relative to the expected separation, and the amplitude was fixed to the theoretical value of [O III]$\lambda$5007/2.98 (Osterbrock, 1989). For the H$\alpha$, $[$N II$]\lambda\lambda$6548,6584 complex, a triple Gaussian was used; the amplitude of [N II]$\lambda$6548 was fixed to be [N II]$\lambda$6584/3.071 (Osterbrock, 1989), the central wavelengths of [N II] were fixed relative to H$\alpha$, and the width of [N II]$\lambda$6548 was fixed to be that of $\lambda$6584. For Seyfert galaxies where a broad H$\alpha$ component was detected (FWHM $\gtrsim$ 1000 km/s), an additional underlying Gaussian was added.

To avoid contamination from AGN in our SFR calculations, and to study the impact that the presence of an AGN has on PAH dust luminosity, we classify our source spectra. We identify AGN candidates using a combination of spectroscopic diagnostics, mid-infrared colors, and X-ray cross-matching. An object is considered to be an AGN if it meets at least one of the four following criteria:

1. 1.

   the spectrum has broad permitted emission lines (i.e., FWHM $\gtrsim$ 1000 km/s) consistent with a Type 1 Seyfert galaxy;

2. 2.

   the object is a Type 2 Seyfert galaxy or LINER identified with the “BPT” line ratio diagram such that log([OIII]/H$\beta$)$>$0.61/(log([NII]/H$\alpha$)-0.02-0.1833$\times$$z$) + 1.2+0.03$\times$$z$ (Kewley et al. (2013), Equation 1);

3. 3.

   for MIR-selected AGN, the AKARI/IRC colors are such that N2-N4$>$0 and S7-S11$>$0 [AB magnitude] (Lee et al., 2007; Shim et al., 2013); and/or

4. 4.

   the object has a Chandra X-ray counterpart within 2.5^′′ (Krumpe et al., 2015).

The left panel of Figure 4 shows the usual BPT-[NII] emission line ratio diagram for our sources. Based on this diagram, there are 215 star-forming galaxies (“BPT-SF”) and 154 AGN candidates (“BPT-AGN”). Based on the MIR colors, there are 779 star-forming galaxies (“MIR-SF”) and 157 AGN candidates (“MIR-AGN”). BPT-classified SF galaxies with N2, N4, S7, and S11 filter detections are 99% consistent with their mid-infrared color classifications; 160 BPT-SF galaxies are also MIR-SF galaxies, while only two BPT-SF galaxies are classified as MIR-AGN. However, the likelihood that a BPT-AGN is classified as a MIR-SF galaxy is high: 121 BPT-AGN are classified as MIR-SF galaxies and 8 are classified as MIR-AGN. This is probably because MIR color selection is stricter, tending to exclude objects with low AGN fraction, and many Seyfert 2’s which the BPT diagram reveals. Thus, an MIR-selected AGN is powerful enough to be certainly identified as a BPT-AGN (left panel). The converse is not always true; some BPT-classified AGN have active nuclei which are too faint to contribute much mid-IR continuum. Both panels are color-coded by their AGN fraction as determined through SED fitting (Section 3.1).

We matched our H$\alpha$ detections with their multi-wavelength photometric counterparts, using a matching radius of the full width at half maximum of the telescope point spread function. When available, we merged photometry from the following telescopes: Chandra (Krumpe et al., 2015), GALEX FUV and NUV (Burgarella et al., 2019), CFHT/MegaCam u^∗, g’, r’, i’, z’ (Oi et al., 2014; Hwang et al., 2007), Subaru/Hyper-Suprime Cam g/r/i/z/Y (Oi et al., 2021), CFHT/WIRCAM J and K_S (Oi et al., 2014), KPNO/FLAMINGOS J and H (Jeon et al., 2014), AKARI/IRC N2, N3, N4, S7, S9W, S11, L15, L18W, L24 (Kim et al., 2012; Murata et al., 2013), Spitzer/IRAC1 (3.6 $\mu$m) and IRAC2 (4.5 $\mu$m) (Nayyeri et al., 2018), WISE 12 $\mu$m and WISE 22 $\mu$m (Wright et al., 2010; Jarrett et al., 2011; Cutri et al., 2021), Herschel/PACS green (100 $\mu$m) and PACS red (160 $\mu$m) (Pearson et al., 2019), and Herschel/SPIRE 250 $\mu$m, 350 $\mu$m, 500 $\mu$m bands (Pearson et al., 2017). We note that some of these datasets do not cover the entire AKARI NEP-Wide field, but our multi-wavelength coverage of the NEP-Deep field is reasonably complete.

We modeled spectral energy distributions (SEDs) with the CIGALE software (Code Investigating GALaxy Emission), which is based on the assumption of energy balance between stellar radiation absorbed by dust in the UV/optical and dust emission output in the mid to far-infrared (Boquien et al., 2019; Yang et al., 2020). This implicitly requires that our view of the extinction is similar to what observers along most other lines-of-sight to the galaxy would also measure. CIGALE uses Bayesian analysis that combines stellar population models, AGN and ISM templates²²2We note that CIGALE includes the option to include user-input nebular emission line fluxes (e.g., H$\alpha$, [N II]), which would further constrain the metallicity and $q_{PAH}$ parameters. However, we do not include emission line fluxes in this work, but recommend this option for future studies., and various libraries to determine the best parameters across all possible spectra. We assumed an exponentially decaying star formation history, Chabrier (2003) IMF, Calzetti et al. (2000) dust attenuation model, Draine et al. (2014) dust emission template, and Fritz et al. (2006) AGN emission model. The dust emission model characterizes PAH emission strength based on the $q_{PAH}$ parameter, which is defined as the fraction of the total dust mass in PAH grains with fewer than 10³ carbon atoms. For reference, average Milky Way dust has $q_{PAH}\approx$ 4.6% (Draine & Li, 2007). We note that the Chandra X-ray data were only used to identify X-ray AGN and were not included in the CIGALE fits.

We define the AGN fraction, $frac_{AGN}$, as the ratio of the AGN luminosity to the total luminosity integrated between 5–10 $\mu$m. For star-forming galaxies, we assume $frac_{AGN}=0$. Otherwise, for AGN candidates and unclassified objects (i.e., objects with incomplete spectroscopic/photometric data), we allowed $frac_{AGN}$ to vary. Figure 4 is color-coded by the AGN fraction. The median AGN fraction for BPT-AGN candidates is 0.05 with interquartile range of 0.01–0.13, compared to a median of 0.33 for MIR-AGN candidates with interquartile range of 0.12–0.48. This difference demonstrates the efficacy of the BPT diagram to select AGN with lower $frac_{AGN}$ than the AKARI color-color diagram.

We used a 32-core processor through Amazon Web Services (AWS) cloud computing³³3https://aws.amazon.com/ to run our models at a rate of $\sim$30,000 models/second. For the final sample, we limit objects to those with $0.5<\text{reduced}\chi^{2}<10$. Objects with reduced $\chi^{2}<0.5$ were excluded in order to avoid models with overfitting. Table 3 summarizes our input SED parameters. In Figure 5, we illustrate the multi-wavelength coverage of our photometric data (purple circles) and show some diverse examples of CIGALE best-fit SEDs which have a range of PAH strengths in the mid-infrared. The light blue region from 2-24 $\mu$m highlights the observed wavelength range of the AKARI/IRC, and the dotted vertical line indicates the observed wavelength for the PAH 7.7 $\mu$m blend. The model spectrum (solid black line) represents the sum over the attenuated stellar emission, dust emission, and AGN emission components. Although we include the contribution from nebular emission when modeling SEDs, we omit it in the figure for simplicity.

The total infrared luminosity, $L(TIR)$, has been extensively used in the literature as a SFR indicator, and is known to correlate with PAH luminosity (Takagi et al., 2010). We calculate $L(TIR)$ for individual objects by shifting the best-fit SED to the rest frame and integrating from 8–1000 $\mu$m. In the Appendix, we test the dependence of FIR photometry on the integrated $L(TIR)$ and find that $L(TIR)$ measurements with and without FIR data are consistent with a root-mean-square (RMS) dispersion of 0.096 dex (Figure 25). Therefore, we estimate that our $L(TIR)$ measurements are accurate to within $\sim$20%.

Our final PAH-derived SFR calibration sample consists of main-sequence and starburst galaxies with H$\alpha$ and/or [O II] emission lines and metallicities measured using the N2, O3N2, and/or O32 line ratio indices. The sample includes 319 galaxies with H$\alpha$ detections and 332 galaxies with [O II] detections. There are a total of 443 unique sources. We summarize the main properties of this sample in Figure 6. The redshift, stellar mass, and total IR luminosity range from 0.048$<z<$1.025, 10^8.7 M_⊙$<$ M_∗ $<$10^11.2 M_⊙, and $10^{8.7}$ L_⊙$<$ L(TIR) $<$10¹² L_⊙, respectively. Table 4 lists the emission line fluxes and physical properties for a portion of galaxies within the calibration sample. The full table is published online.

Figure 7 shows the specific star formation rate (sSFR) as a function of redshift of all galaxy candidates with H$\alpha$ and/or [O II] detections. The sSFR, the mass-normalized star formation rate, is given by Kennicutt (1998) as

where $L(TIR)$ is the total IR luminosity of the best-fit SED integrated from rest-frame 8–1000 $\mu$m and (M_∗/M_⊙)/0.61 is the stellar mass converted from Chabrier (2003) to Kennicutt (1983) IMF. “Main-sequence” galaxies are defined as objects with $13\times t^{-2.2}_{cosmic}\leq$ sSFR [1/Gyr] $\leq 52\times t^{-2.2}_{cosmic}$, where $t_{cosmic}$ is the cosmic time elapsed since the Big Bang in Gyr (Elbaz et al., 2011). “Starburst” galaxies and “quenched” galaxies lie above and below this region, respectively. To classify starburst galaxies, we adopt the $R_{SB}$ starburst parameter from Elbaz et al. (2011), which describes the “starburst excess” relative to main-sequence star formation and is defined as the ratio of the sSFR to the main-sequence sSFR. With their definition, starburst galaxies have $R_{SB}\gtrsim 2$ (blue star symbols), main-sequence galaxies have $0.5\lesssim R_{SB}\lesssim 2$ (black triangle symbols), and quenched galaxies have $R_{SB}\lesssim 0.5$ (red circle symbols).

In Figure 8, we plot the observed N3 and L18W magnitudes as a function of redshift for galaxies within our calibration sample (red circle symbols). Of the 443 star-forming galaxies in the calibration sample, 424 galaxies have detected N3 magnitudes (left panel), and 369 galaxies have detected L18W magnitudes (right). For reference, photometric detections of galaxies in the AKARI/IRC NEP-Deep field with photometric redshifts are shown as black dots (Oi et al., 2014). The blue dashed lines in the figure indicate that our sample is flux-limited down to $\simeq$ 20 AB mag at 3 $\mu$m and $\simeq$ 19 AB mag at 18 $\mu$m up to $z\sim 0.7$, which corresponds to a stellar mass limit of $M_{*}\gtrsim 10^{9.5}$ $M_{\odot}$.

We calculate dust-corrected H$\alpha$ luminosities using three methods: the 24 $\mu$m rest-frame monochromatic luminosity ($\lambda$L_λ(24 $\mu$m)) correction, the Balmer decrement (H$\alpha$/H$\beta$), and the color excess E(B-V) derived with CIGALE. The empirical 24 $\mu$m correction to H$\alpha$ defined by Kennicutt (1998) is given as:

where $L(H\alpha)_{obs}$ is the observed H$\alpha$ luminosity and luminosities are in erg/s. We calculate $\lambda$L_λ(24 $\mu$m) as the rest-frame luminosity density of the best-fit SED model at 24 $\mu$m multiplied by the effective wavelength. The intrinsic H$\alpha$ luminosity is given by:

where $A(H\alpha)$ is the attenuation in H$\alpha$ in magnitudes. For star-forming galaxies with H$\alpha$ and H$\beta$ flux measurements, $A(H\alpha)$ is given by the Balmer decrement (BD):

where (F(H$\alpha$)/F(H$\beta$))_obs is the observed flux ratio, and we assume the Calzetti et al. (2000) attenuation law and an intrinsic flux ratio of 2.86 (Case B recombination and T=10⁴ K; Osterbrock (1989)). The SED-derived $A(H\alpha)$ is given by:

where E(B-V)_lines is the color excess of the nebular lines’ light, and we assume a standard Milky Way extinction curve.

Figure 9 shows a comparison of intrinsic H$\alpha$ luminosities for the galaxies in the calibration sample with H$\alpha$ and H$\beta$ flux measurements. Heavily dust-reddened objects with $A(H\alpha)>3$ mag as calculated via the Balmer decrement are shown by red squares. In these objects, the dust-corrected H$\alpha$ luminosity deviates significantly from those derived via the 24 $\mu$m luminosity and predicted CIGALE extinction, as shown in panels (a) and (b). However, as shown by panel (c), the corrections given by the 24 $\mu$m and CIGALE methods produce similar intrinsic H$\alpha$ luminosities, which is most likely due to their shared dependence on MIR photometry. The scatter between the Balmer decrement correction and 24 $\mu$m luminosity correction has been observed in other studies (e.g., Boselli et al., 2015, Shipley et al., 2016). They attribute the discrepancy as being due to uncertainty in the underlying stellar absorption correction for objects with weak H$\beta$ emission.

To correct the [O II] luminosity for dust extinction, we adopt the calibration given in Kennicutt (1998):

where the luminosities are in erg/s. The [O II] SFR is then given by Kennicutt (1998) as:

assuming a Salpeter IMF. We discuss limitations of using the [O II] luminosity to measure SFR in the Appendix.

We input the rest-frame, best-fit SEDs from CIGALE into PAHFIT (Smith et al., 2007; Smith & Draine, 2012) in order to model the dust emission and continuum from 5-35 $\mu$m, and to extract the PAH luminosities. The input best-fit SEDs include the attenuated stellar emission, nebular emission, and dust emission (i.e., modeled AGN emission is removed for AGN candidates). Then, we resample the model spectrum into the Spitzer/IRS wavelength grid. PAHFIT uses a physical decomposition model to simulate the MIR continuum with a sum of modified blackbodies and starlight, and the PAH dust blends with Drude profiles. Because it does not take AGN emission into account, we only consider SED components from the total attenuated stellar light, nebular emission, and dust emission. For the PAH features, we fixed the central wavelengths and FWHM to their standard values and allowed the amplitudes to vary. Figure 10 shows a typical PAHFIT result for a star-forming galaxy. The integrated intensity of the Drude profile is equal to

where $\lambda_{r}$ is the central wavelength of the PAH feature, $\gamma_{r}$ is the fractional FWHM, and $b_{r}$ is the central amplitude (Smith et al., 2007). For example, the luminosities for the PAH 6.2 $\mu$m and 7.7 $\mu$m emission features are given by:

with amplitudes $b_{r}$ in units of W/m²/Hz, luminosity distance $d_{L}$ in meters, and $c$ in $\mu$m/s.

To determine how well PAH luminosities predict SFR, we correlate the PAH 6.2 $\mu$m and 7.7 $\mu$m luminosities with the reddening-corrected H$\alpha$ and [O II] luminosities. Figure 13 shows the PAH 6.2 $\mu$m and 7.7 $\mu$m luminosity as a function of observed (reddened) H$\alpha$ and [O II] luminosity. The four panels show a positive correlation between PAH luminosity and observed H$\alpha$ and [O II] luminosity, with considerable scatter due to reddening effects which have not yet been taken into account. For reference, the symbols are colored by the logarithm of the starburst intensity, $R_{SB}$. Despite the large scatter, there is a visible trend such that galaxies with higher $R_{SB}$ generally have depreciated PAH luminosity for a given observed H$\alpha$ or [O II] luminosity (or enhanced H$\alpha$ or [O II] for a given PAH luminosity).

Figures 14 shows the PAH 7.7 $\mu$m luminosity as a function of the de-reddened H$\alpha$ luminosity. We fit linear relations of the form log $L(PAH\ 7.7\ \mu m)$ = $c_{0}$ + $c_{1}$ $\times$log $L(H\alpha)$ to the star-forming sample. Quenched galaxies (i.e., $R_{SB}<0.5$; see Section 3.2) are shown for reference as red square symbols but are excluded from the fits. The three panels in Figure 14 show results for different methods of correcting the observed H$\alpha$ luminosity as discussed in Section 3.3. Compared to Figure 13, the scatter between $L(PAH\ 7.7\ \mu m)$ and $L(H\alpha)_{24\ \mu m}$ is significantly reduced; the RMS dispersion of the fit is 0.31 dex. The dotted line in panel (a) represents the linear relation from Shipley et al. (2016) derived from Spitzer/IRS galaxies with solar metallicity at $z<0.4$. For the majority of galaxies, $L(PAH\ 7.7\ \mu m)$ is linearly correlated with $L(H\alpha)_{24\ \mu m}$ (i.e., slope $\sim$ 1). However, we find that a small minority, 47/319 (15%), of the galaxies deviate from the fit by -0.4 dex, contributing to the uncertainty in the normalization and non-linear slope. This deviation is consistent with Calzetti et al. (2007), who attribute it to low metallicity and decline to provide a calibration for 8 $\mu$m emission as a local star formation rate indicator on this basis. Panels (b) and (c) show that the correlations between $L(PAH\ 7.7\ \mu m)$ and $L(H\alpha)$ are non-linear with the Balmer decrement and SED fitting corrections, with the Balmer decrement method resulting in the most scatter.

The majority of the quenched galaxies shown in panel (a) of Figure 14 (red square symbols) lie along or slightly above the linear correlation between $L(PAH\ 7.7\ \mu m)$ and dust-corrected $L(H\alpha)$, which is consistent with the physical scenario where intense starbursts lead to the destruction of PAH dust grains due to UV dissociation (Draine et al., 2007); in an environment with low $R_{SB}$, there are significantly fewer young, massive O and B stars compared to actively star-forming galaxies. Based on a sample of local SDSS galaxies, Smercina et al. (2018) found that the dominant PAH emission in post-starburst galaxies can overestimate the SFR relative to traditional ionized gas tracers in the mid-infrared that are typically present in H II regions (e.g., [Ne II]+[Ne III]). French et al. (2023) similarly found that the TIR luminosity overestimates the SFR relative to H$\alpha$, [Ne II]+[Ne III], and other ionized gas tracers when the star formation history is abruptly truncated, which is believed to be the scenario in post-starburst galaxies.

Figure 15 shows the correlations between the PAH 6.2 $\mu$m and 7.7 $\mu$m luminosities with the dust-corrected H$\alpha$ luminosity for the star-forming sample, with main-sequence and starburst galaxies fit separately. For main sequence galaxies, the linear and unity slope fits for the PAH 7.7 $\mu$m luminosities are consistent with the relations found by Shipley et al. The linear relations for $L(PAH\ 6.2\ \mu m)$ in main-sequence galaxies are also consistent within $\sim$0.25 dex. The figure demonstrates that starburst galaxies have lower PAH luminosities for a given intrinsic H$\alpha$ luminosity. If we assume unity slope in the linear fits, the trend for starburst galaxies has a systematic offset of $\sim$0.35 dex relative to that of main-sequence galaxies for both PAH 6.2 $\mu$m and 7.7 $\mu$m features, suggesting that starburst galaxies are either PAH-deficient for fixed H$\alpha$ luminosity or H$\alpha$-enhanced for fixed PAH luminosity.

The offset between starburst galaxies and main-sequence galaxies is also present in the correlation between PAH luminosity and dust-corrected [O II] luminosity, as shown in Figure 16. Starburst galaxies are systematically lower by a factor of 0.3 dex. We also find that linear fits to the main-sequence sample of the form log $L(PAH)=c_{0}+c_{1}\times$ log $L([OII])$ result in a shallower slope compared to fits with respect to $L(H\alpha)$ such that $L(PAH)\propto L([OII])^{0.9}$.

We further examine the dependence of PAH luminosity on starburst intensity and AGN strength by plotting the log ratio of PAH luminosity to intrinsic H$\alpha$ and [O II] luminosity as a function of PAH equivalent width (EW) in Figure 17. The equivalent width is calculated as the flux of the PAH emission feature (given by PAHFIT) divided by the flux density of the continuum components, which includes the starlight, evaluated at the central wavelength of the feature. We note that it is difficult to constrain the values for $EW(PAH\ 7.7\ \mu m)$ for individual galaxies because the calculation relies heavily on assumptions for the underlying continuum shape and blend of Drude profiles. However, the exact equivalent widths do not affect our overall interpretations. The top row of the figure shows the correlations for star-forming galaxies, color-coded by starburst parameter. The bottom row shows the correlations for AGN candidates, color-coded by AGN fraction. Objects with low PAH equivalent width (i.e., $EW(PAH\ 6.2\ \mu m)\lesssim 1\ \mu$m or $EW(PAH\ 7.7\ \mu m)\lesssim$ $4\ \mu$m) scatter downwards from the $L(PAH\ 6.2/7.7\ \mu m)$ and dust-corrected $L(H\alpha)$/$L([O\ II])$ fits. The figure demonstrates that galaxies with higher starburst intensity ($R_{SB}$) have lower $L(PAH\ 6.2\ \mu m)$ and $L(PAH\ 7.7\ \mu m)$ for a given $L(H\alpha)_{24\ \mu m}$ and that AGN with higher AGN fraction have lower $L(PAH\ 6.2\ \mu m)$ and $L(PAH\ 7.7\ \mu m)$ for a given $L([OII])_{24\ \mu m}$. The offset such that galaxies with high $frac_{AGN}$ have low log $L(PAH)/L([O\ II])$ across all equivalent widths is likely due to the destruction of PAH dust grains by hard radiation from AGN. Given that ionization from AGN limits the efficacy of using the intrinsic H$\alpha$ and [O II] luminosities to trace star formation, we exclude AGN from all SFR calibrations.

To further illustrate the effect of starburstiness on PAH intensity, we show the kernel density estimate and marginal distributions of the PAH 7.7 $\mu$m luminosity to intrinsic H$\alpha$ luminosity as a function of $R_{SB}$ on the left column of Figure 18. The log $L(PAH\ 7.7\ \mu m)$/$L(H\alpha)_{24\ \mu m}$ ratio is anti-correlated with $R_{SB}$. In addition, galaxies with higher $R_{SB}$ are associated with lower PAH 7.7 $\mu$m equivalent width ($\lesssim 4$ $\mu$m), as shown by the orange contours. The Spearman rank-order correlation coefficient between $EW(PAH\ 7.7\ \mu m)$ and $R_{SB}$ is -0.48 with p-value of $2\times 10^{-26}$, indicating a moderate monotonically decreasing correlation. The anti-correlation between log $L(PAH)/L(H\alpha)$ and $R_{SB}$ is consistent with the scenario in which PAH dust grains are destroyed by UV dissociation in compact star-forming regions (Peeters et al., 2004; Murata et al., 2014). For example, Murata et al. (2014) observe a deficit in $\nu L(8)/\nu L(4.5)$ for star-forming galaxies with log $R_{SB}>0.5$, which they interpret as a deficit of PAH emission in starburst galaxies.

The right column of Figure 18 illustrates the correlation between log $L(PAH\ 7.7\ \mu m)$/$L(H\alpha)_{24\ \mu m}$ and $q_{PAH}$ parameter derived from best-fit SED models. The $q_{PAH}$ parameter is highly correlated with PAH equivalent width; the Spearman rank-order correlation coefficient between $EW(PAH\ 7.7\ \mu m)$ and $q_{PAH}$ is 0.8 with a p-value of $1.6\times 10^{-98}$. We note, however, that the strong linear correlation between log $L(PAH\ 7.7\ \mu m)$/$L(H\alpha)_{24\ \mu m}$ and $q_{PAH}$ was to be expected given the correlation between $L(PAH)$ and $q_{PAH}$ based on assumptions in the SED modeling.

In Figure 19, we investigate the dependence of PAH intensity on metallicity. Here, we calculate the PAH intensity as the ratio of $L(PAH\ 7.7\ \mu m)$ to $SFR(H\alpha)_{24\ \mu m}$, both averaged over stellar mass. The stellar masses range from $M/M_{\odot}<10^{9}$ to $M/M_{\odot}>10^{10.5}$ in bins of 0.5 dex. The figure shows the correlation between PAH intensity and average metallicity, $<$12+log(O/H)$>$, calculated using either the N2 (top panel) or O3N2 (bottom panel) index. Main-sequence galaxies (i.e., $0.5<R_{SB}<2$) and starburst galaxies (i.e., $R_{SB}>2$) are shown as circle and star symbols, respectively. Table 5 lists the number of main-sequence and starburst galaxies included in each bin. There is a positive correlation between PAH 7.7 $\mu$m intensity and metallicity in both main-sequence and starburst galaxies, with starburst galaxies having systematically lower PAH 7.7 $\mu$m intensity per stellar mass.

There has been some concern that these “strong-line” ratios might give systematically incorrect metallicities. One possibility, for example, is that the N/O ratio in strongly star-forming galaxies might deviate from the solar value (Masters et al., 2014; Henry et al., 2021; Spinoglio et al., 2022). However, we do not find any clear systematic differences between the metallicities estimated with or without the [OIII] line in Figure 19. Another concern is that electron-temperature based metallicities (using the [OIII]4363 emission line which is too weak to measure in our spectra) might yield systematically lower metallicities than our strong-line methods (Shin et al., 2021; Ly et al., 2016a, b). However, even if a systematic offset of a few tenths of a dex were applied to all of our [O/H] estimates, our finding of a relative increase of PAH strength with metallicity is not significantly changed, other than in the normalization.

In Figure 20, we compare the PAH contribution to the total IR luminosity in star-forming galaxies and AGN by showing the distributions of the total PAH luminosity, $L(PAH)$, to $L(TIR)$ for objects with at least one Herschel FIR detection. We define $L(PAH)$ as the sum of luminosities from PAH emission features at 6.2, 7.7, 8.6, 11.3, 12.6, and 17 $\mu$m. Star-forming galaxies are shown in the top panel and AGN are shown in the bottom panels, with the bottom-most panel separating out “weak” AGN (i.e., $frac_{AGN}<0.1$) from stronger AGN (i.e., $frac_{AGN}>0.1$). We find that $L(PAH\ 7.7\ \mu m)$ contributes a median $\sim 45$% of the total PAH luminosity, and that the total PAH luminosity in turn can contribute up to $\sim 20$% of $L(TIR)$. The median $L(PAH)/L(TIR)$ for star-forming galaxies (main-sequence and starburst galaxies) is 0.08 with an interquartile range of 0.05 to 0.14. The median in the AGN sample is lower at 0.06 with an interquartile range of 0.03 to 0.10. Objects that are more AGN-dominant with $frac_{AGN}>0.1$ have lower $L(PAH)/L(TIR)$ with a median 0.04. For star-forming galaxies, the 10-90% range in $L(PAH)/L(TIR)$ is 0.03 to 0.18, and is plotted as the pale blue shaded region in all three panels in Figure 20. Our results are consistent with those of Shipley et al. (2013). A simple interpretation is that up to half of the TIR luminosity in AGN may be produced by the AGN itself, rather than from star formation. Previous detailed studies of the far-IR suggest that most of this AGN luminosity emerges at wavelengths shortward of 100 $\mu$m (Spinoglio et al., 2002).

In this section, we develop our PAH SFR calibration based on the dust-corrected H$\alpha$ and [O II] luminosities. Given that the AGN fraction is a significant driver of $L(PAH)/L(H\alpha)$ and $L(PAH)/L([OII])$, as demonstrated in Figure 17, our calibration sample excludes AGN and is limited to the 443 unique star-formation-dominated galaxies (Section 3.2). Therefore, $frac_{AGN}$ is not included in the following fits.

Assuming that PAH luminosity corrections for metallicity and starburst intensity are both simply linear, we perform a multi-linear regression to the PAH luminosity with the dust-corrected H$\alpha$ luminosity, R_SB parameter, and metallicity as the independent variables in the form of: log L(PAH) [erg/s] = $c_{0}$ + $c_{1}\times$ log R_SB + $c_{2}\times$(12+log(O/H)_index - 8.6) + $c_{3}\times$[log $L(H\alpha)_{24\ \mu m}$ - 42] [erg/s], where the constant 8.6 is the median metallicity and 42 is the median logarithm of the intrinsic H$\alpha$ luminosity. The relation for predicting the 7.7 $\mu$m PAH luminosity as the dependent variable is: