The Age Dependence of Mid-infrared Emission around Young Star Clusters

The LEGUS sample consists of 50 galaxies between 3.5 and ∼16 Mpc distance. Within this distance range, the Spitzer/IRAC camera, with a ∼2'' point-spread function (PSF) subtends regions between 30 and 150 pc. In order to maximize the number of isolated star clusters within each IRAC PSF, we choose a maximum distance of 5 Mpc for the most distant galaxy we analyze, which subtends regions ∼50 pc in size or the typical size of H ii regions. In addition, we require the LEGUS catalogs to include at least 100 identified star cluster candidates (class 1, 2, or 3; see below) to ensure sufficient statistics. This reduces our sample to the 5 galaxies listed in Table 1. In this table we list, in addition to cluster numbers, also several basic physical properties for each galaxy, such as morphologies, adopted distances, oxygen abundances, and stellar masses. In Figure 1 we display the Galaxy Evolution Explorer (GALEX) 2-color (far-UV and near-UV) images for these galaxies, which are obtained from the GALEX GR6/7 Data Release.¹²

Zoom In Zoom Out Reset image size

Table 1.  Properties of Selected Galaxies

Properties	NGC 7793	NGC 4395	NGC 4449	NGC 1313	NGC 3738
Morph.a	SAd	SAm	IBm	SBd	Im
Dist. (Mpc)b	3.44	4.30	4.31	4.39	4.90
12 + log(O/H)c	8.50	8.19	8.26	8.21	8.10
M* (M⊙)d	$3.2\times {10}^{9}$	6.0 × 108	1.1 × 109	2.6 × 109	2.4 × 108
N(YSC)e	378	108	284	674	164
N(IRAC4 clumps)f	48	17	25	47	9
N(final)g	39	17	16	23	2

Notes.

^aMorphological type as listed in the NASA/IPAC Extragalactic Database. ^bRedshift-independent or flow-corrected redshift-dependent distance in Mpc, taken from Calzetti et al. (2015). ^cCentral oxygen abundances of the galaxies derived from the empirical calibrations (Pilyugin & Thuan 2005; Pilyugin et al. 2012, 2013) or the direct method (Pilyugin et al. 2010, 2012). The values for NGC 7793, NGC 4395, and NGC 1313 are taken from Pilyugin et al. (2014), while the ones for NGC 4449 and NGC 3738 are taken from Pilyugin et al. (2015). ^dStellar masses in M_⊙, taken from Calzetti et al. (2015). ^eNumber of identified YSCs; see Section 2.2. ^fNumber of IRAC4 sources identified as YSC counterparts; see Section 3.2. ^gNumber of IRAC4 sources in the final sample; see Section 4.

Download table as:  ASCIITypeset image

We take the LEGUS catalogs from the public repository at https://archive.stsci.edu/prepds/legus/dataproducts-public.html (Adamo et al. 2017; Cook et al. 2019) as the parent sample, which include photometry in five bands (NUV, U, B, V, I) and physical properties such as stellar mass (M_*), age, and extinction derived from fitting the spectral energy distribution (SED) of each cluster candidate (Calzetti et al. 2015, their Figure 9). The physical parameters are derived using a traditional χ² minimization approach between models and photometric data, as described in Adamo et al. (2010). The more accurate Bayesian approach, in which model parameters are stochastically sampled from cluster evolutionary models, has been applied to the catalogs of only two LEGUS galaxies (Krumholz et al. 2015). For this reason, we use the physical parameters derived from the χ² approach, also called the deterministic approach (Adamo et al. 2017).

The SED fitting for the deterministic catalogs uses two stellar libraries—Padova-AGB (Bressan et al. 1993; Fagotto et al. 1994a, 1994b; Girardi et al. 2000; Vázquez & Leitherer 2005) and Geneva tracks without rotation (Ekström et al. 2012; Georgy et al. 2013; Leitherer et al. 2014)—and three extinction/attenuation laws using photometry from two types of aperture corrections (average based and concentration index (CI) based). The combination of these choices produces 12 catalogs for each LEGUS galaxy (Adamo et al. 2017). Here, we make use of the catalogs constructed from the average aperture corrections and SED fitting utilizing the Geneva tracks and Milky Way extinction law of Cardelli et al. (1989). The choice of aperture correction or model assumption does not have a significant influence on our final conclusions, as briefly discussed in Section 6. Extensive details about the fitting approach and parameter assumptions can be found in Adamo et al. (2017).

The LEGUS catalogs include visually inspected star cluster candidates assigned to "class 1" (symmetric, compact morphology), "class 2" (compact, but elongated), and "class 3" (multipeak). These have been inspected by at least three individuals, after ensuring that each cluster candidate is detected in at least four of the five LEGUS bands with photometric error σ ≤ 0.3 mag. The Adamo et al. (2017) catalogs use a Q value, which is the probability that χ² exceeds a particular value by chance (Press & Teukolsky 2007), to evaluate the quality of the fit; this parameter is larger than 0.1 if the fit is good (Press & Teukolsky 2007; Adamo et al. 2010). Note that the Q parameter is not the reduced χ² (${\chi }_{\mathrm{red}}^{2}$), which is close to 1 for a good fit. To restrict the goodness of fit to the higher quality results, as suggested by Adamo et al. (2010), we only include YSCs with Q > 0.001 (∼80% of cluster candidates with class = 1, 2, or 3) in our sample, indicating that the best-fit models are acceptable (Press & Teukolsky 2007). The resulting sample has a median and 68% scatter of ${0.31}_{-0.28}^{+0.49}$ and ${1.20}_{-0.87}^{+1.89}$ for Q and ${\chi }_{\mathrm{red}}^{2}$, respectively. See Adamo et al. (2017) for more details.

The IRAC1 (3.6 μm) and IRAC4 (8.0 μm) images for our sample galaxies are taken from the Spitzer Local Volume Legacy (LVL; Dale et al. 2009) survey,¹³ which provides IRAC and MIPS images for 258 galaxies out to 11 Mpc. Most of the LEGUS galaxies with public YSC catalogs are also included in the LVL sample by design. The IRAC images are resampled to a pixel scale of 075, and the FWHMs of their PSFs are ∼16 and ∼19 for IRAC1 and IRAC4, respectively (Dale et al. 2009).

The IRAC1 images are convolved to match the PSF of IRAC4 images using the Photutils package (Bradley et al. 2019). Under the assumption that all the 3.6 μm emission is from stars, we remove stellar emission from the IRAC4 images using the formula of Helou et al. (2004). A small fraction (about 5%–33%) of the 3.6 μm flux is due to 3.3 μm PAH emission (Meidt et al. 2012; Robitaille et al. 2012). We ignore this small contribution, which is expected to have negligible impact (≲1.5% for our sample) on the dust-only 8 μm images.

The IRAC4 photometry is performed on the dust-only (star-subtracted) images utilizing the PHOT task in the PyRAF¹⁴ interface of IRAF (Tody 1986, 1993). Due to the large PSF of the IRAC4 images, most of the sources are blended with their neighbors. To extract accurate fluxes for each source, we limit our analysis to relatively isolated sources and perform photometry in small apertures with aperture corrections determined in a way similar to that described in Adamo et al. (2017).

We mark the locations of all the clusters on the cutouts of optical and IRAC4 images and inspect each cutout by eye to identify YSCs with IRAC4 counterparts within 225.¹⁵ For each IRAC4 counterpart, we perform photometry with aperture radii of 3, 3.5, and 4 pixels and sky annuli set at a radius of 10 pixels and a width of 1 pixel. The total fluxes from these three apertures after aperture corrections are in good agreement with each other within the uncertainties. However, we still find that, for some sources, a larger photometric aperture can result in a slightly elevated flux. To avoid any contamination from either neighbors or local background, we adopt 3 pixels (225) as the final aperture radius. Additional examination of the radial profiles of each source is carried out by visual inspection, and IRAC4 sources without well-defined profiles within 3 pixels, due to confusion with neighboring sources, are also removed. In Figure 2, we show two examples of YSCs with IRAC4 counterparts. The top row has only one YSC within the IRAC4 aperture, while the bottom row is an example of multiple clusters within one aperture. The centers of IRAC4 clumps might have small offsets from the coordinates of YSCs in the Adamo et al. (2017) catalogs even in the case of a single cluster within one IRAC4 aperture. Therefore, we allow the PHOT task to recenter before performing photometry. From the radial profile shown in the top row of Figure 2, it is obvious that adopting a photometric aperture of 4 pixel radius could overestimate the 8 μm fluxes due to the contamination of nearby bright sources.

Zoom In Zoom Out Reset image size

Aperture corrections from our 3 pixel aperture radius to an infinite aperture are derived for our sources as follows. We first identify some isolated local peaks on the IRAC4 images that exhibit clean radial profiles and seem not to be significantly contaminated by nearby sources. Taking these selected sources as reference clumps, we then perform photometry on a set of aperture radii (from 0.5 to 10 pixels with a step of 0.5 pixels) for the selected sources to obtain the growth curves. The sky (local background) is measured within an annulus at 10 pixel radius, 1 pixel wide. After visual inspection of all the radial profiles in the reference sample, this aperture is found to be large enough to enclose most of the flux of isolated sources. Therefore, we adopt fluxes within 10 pixel radius (corresponding to ∼145 pc at a distance of 4 Mpc) as total fluxes for the reference clumps that are used to normalize growth curves. An example of the reference clumps, including the IRAC4 image, the radial profile, and the constructed growth curve, in NGC 7793 is shown in the top row of Figure 3. It is evident that this source is isolated enough to obtain a robust growth curve.

Zoom In Zoom Out Reset image size

For each galaxy, at least five reference clumps are used to construct the growth curve, and the mean curve for each is adopted as the final curve. In the bottom row of Figure 3, we present all the growth curves for NGC 7793 (left panel) and all the curves used in this work (right panel). Because the number of isolated local peaks is not sufficiently large to obtain a growth curve for NGC 4449 and NGC 3738, we take the ones from galaxies at similar distances. In other words, NGC 4395 (NGC 1313) shares the same growth curve with NGC 4449 (NGC 3738). The choices of growth curves for NGC 4449 and NGC 3738 only alter the total 8 μm fluxes by up to 20%. Such fluctuations do not significantly change the data and thus the conclusions of this work. For a given aperture, the aperture correction is estimated as the ratio between the fluxes within a 10 pixel radius and given aperture (3 pixels) for the mean growth curve. The standard deviation of the corrections given by all growth curves within one galaxy is also computed and combined with all other uncertainties.

The IRAC Instrument Handbook¹⁶ specifies that sources smaller than 8''–9'', like our sources, should be treated as point sources. The IRAC point-source photometry is calibrated on a 12'' radius aperture, implying that an additional aperture correction needs to be applied to our 75 aperture photometry. Table 4.7 of the Handbook provides a correction of ∼1.05 between a 75 aperture photometry and the default 12'' point-source calibration, which we apply to all our 8 μm measurements.

Given that the 8 μm fluxes are measured within an aperture radius of 3 pixels (225) on the IRAC4 images, all sources enclosed by this aperture may contribute to the 8 μm flux. For this reason, when we consider the total stellar mass corresponding to the measured 8 μm fluxes, the mass contribution from the unidentified clusters or individual stars surrounding the identified YSCs should be included.

We perform photometry on the F814W LEGUS images with the same aperture as the one used for the IRAC4 photometry (i.e., the cyan circles in Figure 2, left) to obtain the total I-band flux in the aperture (${f}_{I,\mathrm{ap}}$). The M_* and I-band fluxes of individual identified YSCs are retrieved from the Adamo et al. (2017) catalogs. We sum the I-band fluxes of all identified YSCs within the aperture to get the total flux counted by our selection (${f}_{I,\mathrm{sel}}$). Assuming that all the objects (star clusters and individual stars) within one IRAC4 aperture have a similar mass-to-light ratio (${M}_{* }/{L}_{I}$), the ${f}_{I,\mathrm{ap}}/{f}_{I,\mathrm{sel}}$ ratio is applied to the total M_* of the identified YSCs within the aperture as a mass correction. The assumption of a constant mass-to-light ratio within each aperture is equivalent to assuming that the bulk of the I-band light originates from same-age stellar populations. This assumption is likely to be acceptable within our small regions (≲50 pc). We show in Appendix B that our regions are better described by more complex stellar populations than a single-age one; however, the impact on the estimated mass in the IRAC4 aperture is small (see Figure B1). For all selected IRAC4 clumps, ${f}_{I,\mathrm{sel}}$ only represents a small fraction of ${f}_{I,\mathrm{ap}}$, leading to a fairly large mass correction with a median and 68% range of ${4.93}_{-2.51}^{+5.89}$.

We first examine how the assumed simple stellar population (SSP) models affect the PAH emission. The ingredients in the SSP models considered here are the stellar evolutionary tracks and metallicities, while the IMF is fixed to the Kroupa (2001) one. We adopt five tracks for the Starburst99 models we generate: GENEVA STD, GENEVA HIGH, PADOVA STD, PADOVA AGB, and GENEVA V00. The first and second tracks are the Geneva models with standard (Schaller et al. 1992; Charbonnel et al. 1993; Schaerer et al. 1993a, 1993b) and high (Meynet et al. 1994) mass-loss rates, respectively. The third and fourth tracks are the original Padova models (Bressan et al. 1993; Fagotto et al. 1994a, 1994b; Girardi et al. 2000) and the ones with thermally pulsing AGB stars (Vázquez & Leitherer 2005), respectively. The last one is the updated Geneva models with zero rotation (Ekström et al. 2012; Georgy et al. 2013; Leitherer et al. 2014). More details can be found in the online document¹⁹ or the relevant literatures (Vázquez & Leitherer 2005; Leitherer et al. 2014). Given that the metallicities of four out of five selected galaxies are about half of the solar metallicity (see Table 1), we only consider three metallicities around this value: Z = 0.004, 0.008, 0.02.

Adopting a dust model with MW3.1_60 (${q}_{\mathrm{PAH}}=4.58 \%$), f_abs = 1, and instantaneous star formation, we run Starburst99 from 0.01 to 500 Myr. f_abs = 1 means that all the optical light is absorbed by dust. This assumption would not be sensible given that clusters in our sample are observed optically. However, we here aim to maximize the PAH emission of each SSP model and thus adopt f_abs = 1 and the highest PAH fraction in the Draine & Li (2007) dust model. We will show later (Sections 5.3 and 5.4) that varying the q_PAH and f_abs values chosen will only change the normalization but not the shape of the relation derived from the models. The predicted relations between $\nu {L}_{\nu ,8}/{M}_{* }$ and the ages for different SSP models are shown in Figure 7. For a fixed SSP model, the PAH emission increases slowly at an ${\mathrm{age}}_{{M}_{* }}$ of ∼1 Myr, reaches its maximum at an ${\mathrm{age}}_{{M}_{* }}$ of ∼3 Myr, and then decreases rapidly toward the high-${\mathrm{age}}_{{M}_{* }}$ end. It can be seen from Figure 7 that changing either the evolutionary track or the metallicity of the SSP models does not alter the normalized PAH emission significantly. Note that there are numerous young sources laying below the (maximal) model predictions plotted in Figure 7, we will discuss this scatter in Sections 5.3 and 5.4.

Zoom In Zoom Out Reset image size

Previous works have demonstrated that there is a positive correlation between the PAH emission/abundance and metallicity (e.g., Engelbracht et al. 2005; Madden et al. 2006; Draine et al. 2007; Smith et al. 2007; Rémy-Ruyer et al. 2015). Although we change the metallicity of our tracks, we do not include the resulting metallicity-dependent PAH abundance, which will be examined in Section 5.3. The conclusion we can draw from Figure 7 is that the change in metallicity in the SSP model does not significantly alter the PAH emission in term of input energy.

Due to the fact that the aperture used for IRAC4 photometry (3 pixels or 225) always encloses more than one stellar object, the observed PAH emission might arise from YSCs with different ages. Although we have applied a criterion that limits the age difference of the YSCs (${\mathrm{age}}_{\max }/{\mathrm{age}}_{\min }\leqslant 3$), it is still possible that L_ν,8 receives additional contributions from young stellar objects that are not in our sample due to their UV-to-I photometry (i.e., detections in less than four bands) or unacceptable SED fitting results (Q ≤ 0.001). Thus, it is reasonable to use complex stellar populations to model the observations, especially for the old objects that lie above all the model predictions plotted in Figure 7.

For simplicity, we consider a continuous star formation model with an SFR of 1 M_⊙ yr⁻¹. It is important to note that the predictions of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation are independent of the SFR assumed here, because ${\mathrm{age}}_{{M}_{* }}$ is weighted by M_* and the PAH emission is normalized by M_*. The Geneva tracks with zero rotation (GENEVA V00 in Starburst99) are adopted as our default tracks. Three metallicities (Z = 0.002, 0.008, 0.014) are used to run the code. Again, a dust model of MW3.1_60 (q_PAH = 4.58%) and a dust absorption fraction of f_abs = 1 are assumed when computing $\nu {L}_{\nu ,8}$ in order to minimize the number of free parameters. In Figure 8, we show the comparison of the predicted PAH emission for different star formation histories (SFHs).

Zoom In Zoom Out Reset image size

As concluded in the previous section, the assumed metallicity in the stellar SED does not have a strong effect on the output total UV-to-NIR luminosity, and thus the PAH emission. For young objects with an ${\mathrm{age}}_{{M}_{* }}\lesssim 5\,\mathrm{Myr}$, the strength of the PAH emission does not change with the SFH. However, for an ${\mathrm{age}}_{{M}_{* }}\gtrsim 5\,\mathrm{Myr}$, a continuous star formation model can greatly enhance the PAH emission, and this enhancement increases with the increasing age of the population. In case of multiple star clusters within one aperture, we find that, at fixed ${\mathrm{age}}_{{M}_{* }}$, the total PAH emission of star clusters with different ages tends to be stronger than that of clusters with similar ages, especially for large ${\mathrm{age}}_{{M}_{* }}$. This trend is also suggested by the sources with ${\mathrm{age}}_{\max }/{\mathrm{age}}_{\min }\gt 3$ (pink hexagons) in the left panel of Figure 4.

In the previous subsections, we only consider the dust model with ${q}_{\mathrm{PAH}}=4.58 \%$ and find that the predicted PAH luminosity only matches a small part of the observations. As shown in Figure 6, varying the q_PAH can significantly change the $\nu {L}_{\nu ,8}/{L}_{\mathrm{TIR}}$ and thus the normalized 8 μm luminosity. Here, we go further to see how the selection of dust models (i.e., ${q}_{\mathrm{PAH}}$) influences the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation.

Assuming our fiducial Geneva tracks with zero rotation and Z = 0.008, we run Starburst99 with both instantaneous and continuous ($\mathrm{SFR}=1\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$) star formation. The model predictions for these two SFHs are shown in Figure 9. The change in q_PAH from 0.47% to 4.58% elevates $\nu {L}_{\nu ,8}/{M}_{* }$ by an order of magnitude as expected. What is not expected is that the range of q_PAH brackets the entire scatter observed in the relation between $\nu {L}_{\nu ,8}/{M}_{* }$ and ${\mathrm{age}}_{{M}_{* }}$. Rémy-Ruyer et al. (2015) indicated that a global q_PAH of 4.58% (the Milky Way PAH abundance) is suitable for these galaxies with subsolar metallicity. Particularly, for NGC 7793, ${q}_{\mathrm{PAH}}={4.30}_{-0.37}^{+0.32}$ was derived in Rémy-Ruyer et al. (2015), while for NGC 4449, the global q_PAH is ∼2%, ${3.02}_{-0.09}^{+0.14}$, and ∼3.2% from Karczewski et al. (2013), Rémy-Ruyer et al. (2015), and Calzetti et al. (2018), respectively.

Zoom In Zoom Out Reset image size

For the instantaneous starburst case, the observations of the youngest populations (${\mathrm{age}}_{{M}_{* }}\lesssim 10$ Myr) can be explained by models with various PAH abundances, whereas for older populations (${\mathrm{age}}_{{M}_{* }}\sim 100$ Myr), some sources exhibit a significant enhancement of PAH emission compared to the model predictions. For the continuous star formation, the predicted $\nu {L}_{\nu ,8}/{M}_{* }$ is in good agreement with the observations up to an ${\mathrm{age}}_{{M}_{* }}$ of ∼200 Myr. These results suggest that additional heating sources of PAH emission enclosed by the IRAC4 aperture might be missed by our YSC selection, especially for sources with large ${\mathrm{age}}_{{M}_{* }}$. Further discussion about sources with excess PAH emission is presented in Section 6.1.

Finally, we fix the SSP model (GENEVA V00 and Z = 0.008), the SFH (instantaneous star formation and continuous star formation with SFR = 1 M_⊙ yr⁻¹), and the dust model (${q}_{\mathrm{PAH}}=4.58 \%$), and vary f_abs from 0.1 to 0.9 (Figure 10). Here, f_abs is the fraction of UV-to-NIR stellar light absorbed by dust and reemitted in the IR, as indicated in Equation (1). Similar to the variations in the dust models, the predicted PAH emission reproduces the observed scatter in the data, while the effect of increasing f_abs is very similar to that of increasing q_PAH.

Zoom In Zoom Out Reset image size

Because we have assumed a fixed U for a given dust model, a constant $\nu {L}_{\nu ,8}/{L}_{\mathrm{TIR}}$ is assigned to each dust model. Therefore, the role of f_abs in determining ${L}_{\nu ,8}$ is the same as that of ${q}_{\mathrm{PAH}}$ although they occur under different conditions. In other words, in the current models, f_abs is completely degenerate with q_PAH. However, given the small range of metallicity variation in our sample of five galaxies, and recalling that the PAH emission is approximately proportional to the metallicity (Rémy-Ruyer et al. 2015), we should expect a small variation in q_PAH. Thus, we tend to favor variations in f_abs as the main reason for the observed scatter in $\nu {L}_{\nu ,8}/{M}_{* }$ at fixed ${\mathrm{age}}_{{M}_{* }}$.

Both Figures 9 and 10 suggest that our observations are better described by an instantaneous burst model for ages below ∼10⁷ yr and by continuous star formation for ages above ∼10⁷ yr, after including the effects of varying ${q}_{\mathrm{PAH}}$ and/or ${f}_{\mathrm{abs}}$. This result is supported by SFH studies based on color–magnitude diagrams. NGC 4449 and NGC 3738 are found to have approximately constant SFHs in the last 200 Myr (Cignoni et al. 2018, 2019; Sacchi et al. 2018). Sacchi et al. (2019) also reported that the SFH of NGC 7793, which contributes 40% of our final sample (see Table 1), is not bursty and suggests inside-out star formation in this galaxy. Further discussion can be found in Appendix B, where we demonstrate that complex stellar populations are still needed to explain the observations of some of our sources.

The Draine & Li (2007) dust model assumes a fixed spectral shape for the local radiation field (i.e., the shape is that of the solar neighborhood, estimated by Mathis et al. 1983) and ignores its variations. As the stellar populations age, the extant model only considers the decrease in the amount of input energy, but ignores the variations in the hardness of the local radiation field. The radiative cross section of the PAHs (e.g., Figure 1 in Draine 2011) is wavelength dependent and increases toward the UV, which means that for YSCs with massive OB stars, and thus more UV photons, the absorption/emission of the PAHs will be enhanced. Because the effect on the absorption is included in f_abs and discussed above, here we only consider the effect on the emission in terms of $\nu {L}_{\nu ,8}/{L}_{\mathrm{TIR}}$, which increases for a harder starlight spectrum.

Draine (2011) reported that $\nu {L}_{\nu ,8}/{L}_{\mathrm{TIR}}$ increases by a factor of 1.57 compared to the Draine & Li (2007) model if the assumed spectrum is replaced by a 2 × 10⁴ K blackbody (cut off at 13.6 eV). This factor will only slightly boost the predicted $\nu {L}_{\nu ,8}/{M}_{* }$ at the small-${\mathrm{age}}_{{M}_{* }}$ end, implying a subtle effect on the overall trend of the predicted ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation. Furthermore, as the starlight travels through H ii regions before heating the PDR, dust in the H ii regions will soften the starlight and reduce the difference between the assumed Mathis et al. (1983) spectrum and the one from a central YSC. Thus, the presence of a harder local radiation field might slightly enhance the PAH emission for YSCs and contribute to the observed scatter, but is not expected to be the main driver of the overall trend observed in Figure 4.

As mentioned in Section 2.2, the YSC catalogs provided by Adamo et al. (2017) include 12 catalogs for each galaxy, resulting from the combination of two aperture-correction methods, two stellar libraries, and three extinction/attenuation schemes. The fiducial ones used in this work are derived from SED fitting with the Geneva tracks without rotation and the Milky Way extinction curve of Cardelli et al. (1989) based on photometry from average aperture correction. Here we test whether the choice of YSC catalogs changes our results.

We take the catalogs generated from the Padova-AGB tracks with the Cardelli et al. (1989) extinction curve and average aperture correction to test the choice of the stellar libraries. Based on the IRAC4 counterparts identified visually, we redo all the calculations except for the IRAC4 photometry. The median and corresponding 1σ scatter of the ratios between the new values and the fiducial ones are ${1.00}_{-0.03}^{+0.02}$ and ${1.00}_{-0.12}^{+0.03}$ for $\mathrm{log}\,{\mathrm{age}}_{{M}_{* }}$ and $\mathrm{log}(\nu {L}_{\nu ,8}/{M}_{* })$, respectively. The overall distribution of these sources on the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ plane remains nearly unchanged. The model predictions of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relations with varying ${q}_{\mathrm{PAH}}$ or ${f}_{\mathrm{abs}}$ still cover all the observations for continuous star formation, while two more sources with ${\mathrm{age}}_{{M}_{* }}\gt 40\,\mathrm{Myr}$ move beyond the predicted curve of ${q}_{\mathrm{PAH}}=4.58 \%$ for instantaneous star formation. Therefore, the choice of evolutionary tracks in the SED fitting has no significant impact on our results.

Similarly, replacing the Cardelli et al. (1989) extinction with the starburst attenuation law of Calzetti et al. (2000) results in small changes around the fiducial values in either ${\mathrm{age}}_{{M}_{* }}$ or $\nu {L}_{\nu ,8}/{M}_{* }$ for the two extinction assumptions (i.e., that the ionized gas has either the same or higher attenuation than the stars).

With regard to the aperture-correction method, the average one changes the normalization of the SED of the YSCs but not the shape, thus only M_* is affected, and age and extinction remain mostly unchanged after the correction (Adamo et al. 2017). Conversely, the CI-based correction changes both the normalization and the shape of the SED, so that M_*, age, and extinction are all affected. However, these physical properties derived from the CI-based aperture correction show overall agreement with those derived from the average-based one, with some scatter (Adamo et al. 2017; Cook et al. 2019). When we use the catalogs with the CI-based aperture correction, a slightly larger scatter around the fiducial values is found for both ${\mathrm{age}}_{{M}_{* }}$ and $\nu {L}_{\nu ,8}/{M}_{* }$. The median and the 1σ scatter are ${1.00}_{-0.05}^{+0.02}$ and ${1.00}_{-0.06}^{+0.11}$ for ${{\rm{log}}{\rm{age}}}_{{M}_{\ast }}$ and $\mathrm{log}(\nu {L}_{\nu ,8}/{M}_{* })$, respectively. Such scatters do not significantly change the overall distribution of our sources on the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ plane. Therefore, our conclusions drawn from the fiducial catalogs are not impacted by the choice of YSC catalogs.

In the last section, we have demonstrated that both ${q}_{\mathrm{PAH}}$ and ${f}_{\mathrm{abs}}$ play similar roles in predicting the PAH emission. In principle, under the frame of instantaneous star formation, the combination of these two parameters can explain all the observations below the model predictions of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation presented in Figure 9, even for the two we labeled as outliers in Section 4. For example, assuming ${q}_{\mathrm{PAH}}=2.50 \%$, ${f}_{\mathrm{abs}}$ should be 0.032 and 0.003 to match the observations of NGC 1313-E 640 and NGC 4395-S 622, respectively.

However, for the IRAC4 sources located above the model predictions, it is difficult to reconcile them with our simple model. In Section 5.3, we argue that such PAH excess indicates additional heating sources that might be missed by our YSC selection. In order to understand this PAH excess, we mark three 8 μm sources located above the model predictions of instantaneous star formation as colored open circles in the top panel of Figure 9. The Hα, optical RGB, and IRAC4 cutouts of these sources are shown in Figure 11. Two of them are from NGC 1313, while the other one is from NGC 4449. The Hα images are retrieved from the HLA: F657N for NGC 1313 (Program: GO–13773; PI: Rupali Chandar) and F658N for NGC 4449 (Program: GO–10585; PI: Alessandra Aloisi).

Zoom In Zoom Out Reset image size

We first attempt to use the identified YSCs alone to explain the PAH emission excess by performing SED fitting with additional Hα photometry. There are four clusters within the IRAC4 apertures plotted in Figure 11. We measure Hα fluxes of these clusters following the photometry method described in Adamo et al. (2017). Due to the comparable spatial resolution of the Hα images to the broadband images of the LEGUS survey, the isolated star cluster samples identified by Adamo et al. (2017) are adopted to construct the growth curves and thus the average aperture corrections. Small offsets between the coordinates of the Hα and LEGUS images recovered by the World Coordinate System information are corrected utilizing IRAF tasks. We have applied our photometry code to the LEGUS images and find good agreement between our measurements and those retrieved from the Adamo et al. (2017) catalogs.

The model photometry for SED fitting is generated from the Yggdrasil population synthesis code (Zackrisson et al. 2011) in which a Kroupa (2001) IMF is assumed. To match the catalogs adopted, we use the Geneva tracks without rotation, together with metallicities of Z = 0.004, 0.008, and 0.02. The Milky Way extinction curve of Cardelli et al. (1989) is adopted, while $E(B-V)$ varies from 0 to 1.5 in steps of 0.01. For nebular emission, a covering factor of 0.5 is assumed. A ${\chi }^{2}$ minimization is implemented to obtain the best-fit model.

We plot the results from the newly fitted M_* and age, which include Hα measurements, in the top panel of Figure 9 as filled circles with error bars with the same colors as the corresponding open circles. Adding the narrowband fluxes in Hα does not significantly change M_* and age for three out of the four clusters. For the first IRAC4 source plotted in Figure 11 (i.e., NGC 1313-W 1871), two identified YSCs are enclosed by the aperture. Comparing with the best-fit results extracted from the Adamo et al. (2017) catalogs, the central one (NGC 1313-W 1871, marked as red circle) has a similar M_* (${8594}_{-1046}^{+1686}\,{M}_{\odot }$ and ${9376}_{-851}^{+1645}\,{M}_{\odot }$) and age (${200}_{-0}^{+100}$ and ${400}_{-200}^{+0}$ Myr), while the one near the aperture edge (NGC 1313-W 1894, marked as blue circle) shows great changes in both M_* (from ${40570}_{-29020}^{+3020}\,{M}_{\odot }$ to ${6356}_{-215}^{+440}$ M_⊙) and age (from ${100}_{-85}^{+0}$ to ${7}_{-0}^{+0}$ Myr). The inclusion of Hα fluxes in the SED fitting results in a slightly larger ${\mathrm{age}}_{{M}_{* }}$ but a smaller total M_* compared to the ones based on the Adamo et al. (2017) catalog. According to the new fitting results, NGC 1313-W 1894 (${\mathrm{age}}_{\mathrm{best}}=7$ Myr) is much younger than NGC 1313-W 1871 (${\mathrm{age}}_{\mathrm{best}}=400$ Myr). Thus, the 8 μm luminosity of that source should be dominated by NGC 1313-W 1894. If we attribute all 8 μm emission to NGC 1313-W 1894, the new location of this source is marked by a blue filled circle and an error bar in the top panel of Figure 9, which is consistent with the model prediction within the uncertainty. Therefore, with the aid of additional narrowband photometry, the PAH excess of NGC 1313-W 1871 can be explained by a much younger YSC within the same IRAC4 aperture.

However, for the other two sources (NGC 1313-E 116 and NGC 4449 885), the new fitting results do not significantly alter the original locations. Even if we use the photometry applying the CI-based aperture correction, such PAH excess still exists. NGC 1313-E 116 has an ${f}_{U,\mathrm{sel}}/{f}_{U,\mathrm{ap}}$ value of 10%, which just passes the second criterion listed in Section 4. If the CI-based photometry is used, ${f}_{U,\mathrm{sel}}/{f}_{U,\mathrm{ap}}$ is only 7.6%, indicating that a large fraction of U-band flux is contributed by other stellar objects within the same IRAC4 aperture, and this source would be removed from our final sample. For NGC 4449 885, the CI-based photometry results in a larger age and a smaller $\nu {L}_{\nu ,8}/{M}_{* }$ value, but it still lies above the model predictions significantly. Therefore, the identified YSCs alone cannot explain the observed PAH excesses. The ages of the two star clusters covered by these two IRAC4 sources are both ∼200–300 Myr, and our measurements should not be as high in $\nu {L}_{\nu ,8}$ as observed. We speculate that additional younger (and, possibly, dust-embedded) clusters missed by our YSC selections or detections in the Adamo et al. (2017) catalogs might account for the PAH excess.

To demonstrate this possibility, we calculate the model predictions of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation for the combination of two YSCs based on the relation of an instantaneous starburst model with ${q}_{\mathrm{PAH}}=4.58 \%$ (Figure 9). The age of the first cluster is set to 300 Myr, which is the same as the two clusters with PAH excess, while the second one has an age that varies from 1 to 500 Myr. The mass ratio of the new clusters changes from 0.01 to 100. In Figure 12, we show how the resulting relations vary with different mass ratios. Clearly, the combination of a 300 Myr old cluster with a younger cluster (i.e., the predictions before the red star) will lead to a significant excess of PAH emission, irrespective of the mass ratio. For instance, in the case of ${M}_{* ,2}/{M}_{* ,1}=0.01$, young individual stars/stellar associations with an age of 5 Myr (red "x" in Figure 12) and ${M}_{* }\approx 60\,{M}_{\odot }$ and 190 M_⊙ are massive enough to boost the PAH emission beyond the model significantly for NGC 1313-E 116 and NGC 4449 885, respectively. Several blue point-like sources can be identified within the IRAC4 apertures of NGC 1313-E 116 and NGC 4449 885 by inspection of the optical RGB images in Figure 11. Therefore, we expect that such blue stellar objects or others with heavy obscuration might contribute to the observed PAH emission excesses of these two sources.

Zoom In Zoom Out Reset image size

Such an excess of PAH emission due to nearby young stellar objects might also happen to YSCs. However, due to the similar age and ${M}_{* }/{L}_{I}$ to the main clusters, the effect would be naturally included by the mass correction. If the missed stellar objects are older than those used in age-dating, the mass correction becomes problematic and thus the PAH emission excess arises. As a result, the PAH emission excess should be more serious toward high-${\mathrm{age}}_{{M}_{* }}$, which is what we observe in Figure 9.

Our 8 μm photometry is performed within the aperture with a radius of 225, which corresponds to a physical radius of ∼30–50 pc for our selected galaxies. For less massive star associations or clusters, we believe that this aperture is large enough to enclose most of the 8 μm fluxes after applying the aperture correction (e.g., Hunt & Hirashita 2009; Lawton et al. 2010). However, for those super star clusters that are powerful enough to ionize a large H ii region, the fixed ∼50 pc aperture might miss the 8 μm flux at some level, leading to an underestimation of $\nu {L}_{\nu ,8}$. For instance, in the most active H ii region in the Large Magellanic Cloud (30 Doradus), the distance between PDR clouds and the main ionizing cluster ranges from ∼11 to 80 pc (Chevance et al. 2016). Given the current data, the fraction of 8 μm flux missed by the fixed aperture is very difficult to estimate and completely degenerate with ${f}_{\mathrm{abs}}$. The effect of missing the 8 μm flux is equivalent to that of varying ${f}_{\mathrm{abs}}$ and might be another possible reason for the observed scatter of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation.

The mass correction described in Section 3.3 assumes that stellar objects (star clusters and field stars) within the same IRAC4 aperture are same-age stellar populations. If the ages of the undetected stellar objects (e.g., individual field stars or undetected star clusters) are significantly different from the detected ones, the mass correction might be problematic.

Here, we consider two cases: (a) the detected star clusters are much younger than undetected stellar objects; (b) the detected star clusters are much older than undetected stellar objects. In the first case, only the detected YSCs contribute to the PAH emission, while the total M_* related to the PAH emission will be overestimated after mass correction due to the existence of field stars. Thus, the undetected old stellar objects might be another reason for the observed scatter of the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation. The two outliers shown in Figure 5 might belong to this case. In the second case, the observed PAH emission excited by the undetected objects might be attributed to the old star clusters, leading to a mismatch of the excited sources of PAH emission and the derived stellar properties. In fact, this should be the reason for the observed PAH emission excess discussed in Section 6.1, but only affects those sources with large ${\mathrm{age}}_{{M}_{* }}$. In Appendix B, we discuss extensively the effects of the undetected stellar populations via the SED fitting using the total fluxes within the aperture.

We have demonstrated that the ${\mathrm{age}}_{{M}_{* }}$–$\nu {L}_{\nu ,8}/{M}_{* }$ relation strongly depends on the PAH abundance and the dust absorption fraction. Even after fixing the PAH abundance using metallicity constraints, it is very difficult to predict the amount of PAH emission for clusters at a given age. It is, however, noteworthy that the model predictions with the highest PAH abundance set an excellent upper limit for the observations, especially for star clusters younger than 10 Myr. Therefore, we fit the model relations derived from a combination of the Geneva tracks with zero rotation and Z = 0.008, together with q_PAH = 4.58% and f_abs = 1, to provide an upper limit to the PAH emission at a given age. Setting $x=\mathrm{log}({\mathrm{age}}_{{M}_{* }}[\mathrm{yr}])$ and $y=\mathrm{log}(\nu {L}_{\nu ,8}/{M}_{* }[{L}_{\odot }/{M}_{\odot }])$, we fit the curves with a piecewise function comprising one constant and two second-order polynomials using the Python version of MPFIT²⁰ (Markwardt 2009). The best-fit results for ${\mathrm{age}}_{{M}_{* }}$ in the range 4–200 Myr are

$$\begin{eqnarray}&&y=\left\{\begin{array}{l}2.0678,\qquad \quad 4\lt x\lt 5.6853,\\ 9.9586-2.8227x+0.2524{x}^{2},\\ \qquad \qquad \quad 5.6853\lt x\lt 6.4786,\\ 19.8870-3.8614x+0.1761{x}^{2},\\ \qquad \qquad \quad 6.4786\lt x\lt 8.3010\end{array}\right.\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&y=\left\{\begin{array}{l}2.0960,\qquad \quad 4\lt x\lt 6,\\ -35.4200+11.9823x-0.9549{x}^{2},\\ \qquad \qquad \quad 6\lt x\lt 6.6839,\\ 5.3867-0.2282x-0.0415{x}^{2},\\ \qquad \qquad \quad 6.6839\lt x\lt 8.3010\end{array}\right.\end{eqnarray} \tag{ 3 }$$

for the instantaneous and continuous star formation models, respectively. These two functions recover the model predictions within 1%.

In Figure 13, we show these best-fit results, overplotted with our sample. The derived function for the continuous star formation model agrees well with the upper limit of our data, while the one from the instantaneous star formation model also gives a good description for the upper boundary of the data. Thus, Equations (2) and (3) can be used to estimate the maximum of PAH emission for either a single star cluster or a star-forming region with either constant or instantaneous SFR at a given age. We stress that these relations are independent of stellar mass and SFR, provided that the measured $\nu {L}_{\nu ,8}$ is normalized by the underlying stellar mass.

Zoom In Zoom Out Reset image size