Evidence for the Accretion of Gas in Star-forming Galaxies: High N/O Abundances in Regions of Anomalously Low Metallicity

Our understanding of the evolution of galaxies is that they continually grow with time. Part of this growth is through mergers between galaxies, but the majority of growth is through a more steady accretion of gas from the cosmic web (e.g., Somerville & Davé 2015). In galaxies currently forming stars, much of this accreted gas ultimately fuels the continuing formation of stars.

Despite the critical role played by accretion, direct observational evidence of significant inflows of gas into galaxies remains scarce outside the Milky Way (Putman 2017). The so-called down-the-barrel observations of intense star-forming galaxies using interstellar absorption lines as probes overwhelmingly reveal gas being expelled from galaxies (Veilleux et al. 2020), with only a small minority of galaxies showing inflows (Martin et al. 2012; Rubin et al. 2012; Zheng et al. 2017). This suggests that either the inflowing gas is not in the warm ionized phase being probed in absorption, or that it covers only a small solid angle as seen from the galaxy.

In this paper we report the results of a complementary approach that provide evidence for on-going accretion of gas based on its chemical abundances. This study is a follow-up to our previous analysis of Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) data (Hwang et al. 2019) in which we found that late-type galaxies in the low-z universe commonly have regions within them where the oxygen abundance (12 + log(O/H)) in the interstellar medium (ISM) is significantly lower than expected (relative to the tight empirical correlation between the oxygen abundance and the local stellar surface mass density for a fixed stellar mass). We defined these as anomalously low-metallicity (hereafter ALM, see Section 2.2) regions, and showed that they correspond to regions of higher than average rates of star formation. We speculated that these ALM regions are sites in which low-metallicity gas has been recently accreted, gotten mixed with the pre-existing more metal-rich gas, and triggered star formation.

The goal of this paper is to test this idea by measuring both O/H and N/O in the ALM regions to search for the tell-tale chemical signature of such an event.

The key to our approach is the behavior of the ratio of nitrogen to oxygen (N/O) as a function of oxygen to hydrogen (O/H). Elements produced by nucleosynthesis are defined as primary or secondary. Primary elements are those whose yield is independent of the initial chemical composition of the star, while the yield of secondary elements depends on the initial composition of the star. Oxygen is a prototypical primary element, produced by helium burning in massive stars (>10M_⊙) and released promptly back to the ISM by core-collapse supernovae. Nitrogen production is more complex (e.g., Arnett 1996; Kobayashi et al. 2020). It is produced from carbon and oxygen as part of the CNO cycle, primarily in intermediate-mass stars (∼4–8 M_⊙), and ejected into the ISM during the post-main-sequence asymptotic giant branch phase. In a low-metallicity star, nitrogen is produced from carbon and oxygen that were previously created via helium burning in that star. In this case nitrogen behaves like a primary element, resulting in a fixed ratio of N/O. However, in higher-metallicity stars, most of the carbon and oxygen used to create nitrogen were pre-existing in the ISM from which the star formed. In this case, nitrogen behaves like a secondary element, with N/O increasing with O/H (e.g., Matteucci & Tosi 1985; Arnett 1996; Henry et al. 2000; Meynet & Maeder 2002). This pattern can be seen in data of chemical abundances of ionized gas in local star-forming galaxies (Garnett 1990; Vila-Costas & Edmunds 1992; Henry et al. 2000; Amorín et al, 2010; Andrews & Martini 2013; Pérez-Montero et al. 2013; Vincenzo et al. 2016; Belfiore et al. 2017). Andrews & Martini (2013) derived a schematic form for N/O versus O/H using SDSS star-forming galaxies (their Figure 14): for 12 + log(O/H) less than 8.5, log(N/O) is constant ∼−1.4; while at higher values of 12 + log(O/H), they found N/O ∝ O/H^1.7.

This locus of N/O versus O/H therefore traces the typical chemical evolution of the interstellar gas polluted by the by-products of stellar nucleosynthesis. In a simple closed-box model, a star-forming region in a galaxy evolves along this locus as more and more gas is converted into stars and the ISM metallicity rises. Eventually, the ISM metallicity becomes high enough for the significant production of secondary nitrogen. There are a number of ways in which a region in a galaxy could depart from this relation between O/H and N/O. For instance, in regions with a significant population of Wolf–Rayet (WR) stars, the gas can be polluted by nitrogen-rich winds from these stars (Walsh & Roy 1989; López-Sánchez et al. 2007; Brinchmann et al. 2008). Fountain-flow models (Bregman 1980; Fraternali & Binney 2008) in which outflowing metal-rich gas from the central region of a galaxy is accreted and mixed with metal-poor gas in the outer region can also produce gas that departs from the normal N/O versus O/H relation (Belfiore et al. 2015).

Of particular interest to the nature of the ALM regions, models in which infalling metal-poor gas is mixed with ambient metal-rich gas (Köppen & Hensler 2005) predict that the mixed gas will be characterized by having an unusually large N/O ratio for a given O/H value. Our specific goal in this paper is therefore to measure the locations of the ALM regions in the N/O versus O/H plane. If, as suggested in Hwang et al. (2019), ALM regions represent sites where accretion of metal-poor gas has occurred, we should see an N/O excess. This technique has recently been used by Loaiza-Agudelo et al. (2020) to show that accretion of metal-poor gas is likely to be responsible for triggering intense star formation in low-z analogs to Lyman break galaxies.

In Section 2 we briefly describe the data used in our study, and we describe our approach to calculate the element abundances in Section 3. In Section 4 we present our results and briefly discuss their implications in Section 5. We summarize our conclusions in Section 6.

Following the same criteria as in Hwang et al. (2019), we focus on late-type galaxies (selected using fracdeV < 0.7, where fracdeV is the fraction of fluxes contributed from the de Vaucouleurs profile), pure star-forming spaxels (selected based on the [N ii] BPT diagram and the Kauffmann et al. 2003 demarcation line) with deprojected local stellar surface mass density > 10⁷M_⊙ kpc⁻², total stellar mass > 10⁹M_⊙, and axis ratio from elliptical Petrosian fitting (b/a) > 0.3 (to exclude edge-on galaxies). We also require signal-to-noise ratio (S/N) > 10 for the emission lines used in BPT diagrams and metallicity calculations.

Hwang et al. (2019) used the strong line calibrator [O ii]3727,3729/[N ii]6584 (O2N2) from Dopita et al. (2013) to estimate the observed (12 + log(O/H))_obs of each spaxel. They computed the expected metallicity for each spaxel (12 + log(O/H))_exp via the interpolation of the tight relationship between local metallicity and the local stellar surface mass density at a given stellar mass (Σ_*–Z relation, Barrera-Ballesteros et al. 2016). They then used the metallicity deviation Δlog(O/H) = (12 + log(O/H))_obs − (12 + log(O/H))_exp to define the ALM spaxels. Here we follow the same approach, but with a different strong line metallicity calibrator: RS32 = [O iii]5007/Hβ + [S ii]6717,6731/Hα from Curti et al. (2020). The photoionization models for the O2N2 indicator assume a relation between the O/H and N/O abundance ratios (Dopita et al. 2016). Since we will compare O/H and N/O of selected ALM spaxels, we need a calibrator that does not involve nitrogen. The RS32 calibrator is a good choice: it has only a mild dependence on the ionization parameter, and is insensitive to dust extinction due to the proximity of the emission line wavelengths.¹⁶ Curti et al. (2020) calibrated the RS32 diagnostic from electron temperature (T_e) based measurements of oxygen abundance in individual galaxies and stacked spectra, and we use the polynomial fitting results in their Table 2 to calculate (12 + log(O/H))_obs here.

Our measured Σ_*–Z relation is shown in Figure 1. Each circle represents the mode of the metallicity distribution inside each stellar mass and surface mass density bin. The step of the stellar mass bin is 0.5 dex, except for the bin of 10¹⁰–10¹¹M_⊙, because there are fewer galaxies with large stellar mass and because the mass-dependence of O/H is weak over this range (Hwang et al. 2019). The bin of M_* = 10^8.5–10⁹ M_⊙ is only used for interpolation, and we do not further analyze the galaxies in this bin. The distribution of the metallicity deviation Δlog(O/H) is shown as the blue histogram in Figure 2. We observe the same low-metallicity tail as in Hwang et al. (2019). To characterize the "normal" spaxels, we fit a Gaussian profile to the positive side of the metallicity deviation distribution, assuming a symmetric negative side. The best fit has σ_Z = 0.037 and is shown as the orange dashed line in Figure 2. The residual (black line in Figure 2) is constructed by subtracting the fitted Gaussian from the negative side of the total distribution. The vertical dotted line marks where the metallicity deviation is three times the scatter. We define ALM spaxels as spaxels having ${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})\lt -3{\sigma }_{Z}$ (−0.111 dex), which makes sure the ALM sample does not include substantial "normal" spaxels characterized by the orange Gaussian profile. The control sample (non-ALM) then consists of spaxels where $| {\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})| \leqslant {\sigma }_{Z}$ (shaded in purple in Figure 2), meaning that they follow the local Σ_*–Z relation well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Extinction corrections for emission line fluxes obtained from DAP assume the average Milky Way extinction curve in Fitzpatrick et al. (2019) Table 3 (with R(55) = 3.02, corresponding to the Milky Way mean value of R(V) = 3.10), and the intrinsic flux ratio of Hα/Hβ = 2.86 (assuming T_e = 10⁴ K, Osterbrock & Ferland 2006). We do not correct for extinction in the rare cases where the measured Hα/Hβ < 2.86. Due to the typically small effect of the Milky Way reddening and the small redshifts of our sample galaxies, this extinction correction should be able to correct for extinction from both the Milky Way and the target galaxy.

We analyzed both the individual spectrum of each spaxel and spectra created by stacking many spaxels in 0.1 dex bins of (12 + log(O/H))_obs. By stacking the spectra we are able to increase the S/N and thus detect the weak auroral lines required to apply the T_e-based direct method to measure O/H and N/O. However, the resulting systematic uncertainties prevent us from using the direct method metallicities to investigate the N/O versus O/H relation (see Appendix A.5). Therefore, we continue to use the strong line diagnostics to determine element abundances.¹⁷ The stacking procedure and direct method calculations are discussed in the Appendix.

For each spaxel, we first calculate the 12 + log(O/H) from the RS32 diagnostic as discussed in Section 2.1. We then calculate N/O using the empirical relation found in Loaiza-Agudelo et al. (2020):

$$\begin{eqnarray}&&\mathrm{log}({\rm{N}}/{\rm{O}})=0.73\times {\rm{N}}2{\rm{O}}2-0.58,\end{eqnarray} \tag{ 1 }$$

where N2O2 = [N ii]6584/[O ii]3727, 3729 (corrected for extinction). Our estimations of O/H and N/O are quite independent as different lines are used. This helps to minimize any potential bias in our comparison of O/H and N/O. Note we are assuming here that the N/O ratio is the same as N⁺/O⁺. This is a safe assumption since the ionization potentials for O⁺ and O⁺⁺ are very similar to those of N⁺ and N⁺⁺ (Garnett 1990).

There exist many strong line metallicity calibrators in the literature, and for the oxygen abundance ratio O/H, it is known that different metallicity calibrators can give results that differ by ∼0.4 dex (Kewley & Ellison 2008). The nitrogen abundance ratio N/O is generally derived from the strong line ratio N2O2. The N/O relation used in the above section is derived based on Lyman break analogs, which are starburst galaxies with high specific star formation rate (sSFR), and our ALM spaxels also have high sSFR (Hwang et al. 2019). Moreover, this N/O calibrator produces very similar results to other N/O relations in recent papers (e.g., Strom et al. 2017).

While we explain why we choose to use RS32 in Section 2.2, we here emphasize that our main result as shown in Figure 3 still holds with combinations of different O/H calibrators and N/O calibrators (see the test figure in Appendix B). The key point here is to be consistent with the calibrator used throughout (from the Σ_*–Z relation calculation and sample selection to exploring the N/O versus O/H relation), since our analysis depends on the relative values of metallicities. We stick to our choice of RS32 for O/H because the R23 and O3S2 calibrators depends strongly on the ionization parameter, and they also show greater scatter around the fitted relation in the range of metallicity where most of our star-forming spaxels lie (Figure 1 in Curti et al. 2020). Our main results are thus robust under different choices of the metallicity calibrators.

Zoom In Zoom Out Reset image size

We plot in Figure 3 the measured N/O versus O/H for the ALM (${\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})\lt -3{\sigma }_{Z}$) and non-ALM ($| {\rm{\Delta }}\mathrm{log}({\rm{O}}/{\rm{H}})| \leqslant {\sigma }_{Z}$) spaxels. The non-ALM spaxels clearly show the trend of N/O increasing with O/H, and the best-fit line is

$$\begin{eqnarray}&&\mathrm{log}({\rm{N}}/{\rm{O}})=2.98\times [12+\mathrm{log}({\rm{O}}/{\rm{H}})]-26.82.\end{eqnarray} \tag{ 2 }$$

The ALM spaxels partially overlap with the non-ALM spaxels, and show a plateau of N/O ∼ −1.4 at low O/H (≲8.5) as expected. We thus take the log(N/O) value from Andrews & Martini (2013) for 12 + log(O/H) ≤ 8.5:

$$\begin{eqnarray}&&\mathrm{log}({\rm{N}}/{\rm{O}})=-1.43.\end{eqnarray} \tag{ 3 }$$

As can be seen, there is a clear offset between the loci of the ALM and non-ALM spaxels. The ALM spaxels lie systematically off the normal relationship, having abnormally large N/O for a given O/H. We interpret this as a signature of the mixing of metal-rich gas (on the linear part of the N/O versus O/H relation) with metal-poor gas (on the flat part of the N/O versus O/H relation). To further quantify this, we calculate log(N/O) and 12 + log(O/H) for a hypothetical gas mixing scenario: we assume that pre-existing gas in the disk with different initial metallicities (following Equation (2)) mixes with accreted low-metallicity gas with half the initial O/H ratio and log(N/O) = −1.43 (following Equation (3)). These abundances in the accreted gas are consistent with measurements of the circumgalactic medium (Prochaska et al. 2017), which represents the likely gas source for accretion. We consider the case where the accreted gas mass is 40% of the total gas mass. The resulting metallicities of mixing are summarized in Table 1 and shown as green arrows in Figure 3. Note that this is meant to be a simple illustrative example of how this mixing process works, and different choices could be made for the chemical composition of the accreted gas.

This mixing model can be further tested by using (12 + log(O/H))_exp values calculated in Section 2.2 to connect the origins of the gas prior to accretion (lying along the non-ALM fit line) to its present location in the N/O versus O/H plot. That is, we assume that the initial O/H value of an ALM spaxel prior to mixing is the value expected based on the Σ_*–Z relation (Figure 1) and that the initial value of N/O is given by Equation (2). For each such initial O/H value, we locate all the ALM spaxels whose (12 + log(O/H))_exp is within 0.02 dex of this initial value, and take the center of their distribution in the N/O versus O/H plot as the ending point of the mixing process. The resulting values for O/H and N/O are summarized in Table 1. In Figure 4, we plot ALM spaxels within the contour that encloses 90% of the total data points, color-coded by (12 + log(O/H))_exp, and show arrows connecting the initial points (same as those in Figure 3) to the ending points. We see that the trajectories of the ALM spaxels are in qualitative agreement with the simple mixing model in Figure 3. This provides further support for the interpretation of the ALM spaxels as regions where pre-existing higher-metallicity gas mixes with accreted lower-metallicity gas.

Zoom In Zoom Out Reset image size

To test the robustness of these results, we subdivide the ALM sample into three different pairs of smaller samples, based on: (1) the equivalent width of the Hα emission line, (2) the stellar mass of the galaxy, and (3) whether the galaxy was classified as merger/close pairs or isolated in Hwang et al. (2019), where it was found that the incidence rate of ALM spaxels is correlated with all these properties. We plot the sub-samples for comparison in Figure 5. Each sub-sample shows enhanced N/O at fixed O/H compared to the expected distribution with no clear differences between the sub-samples of each analyzed pair.

Zoom In Zoom Out Reset image size

One particularly interesting property of the ALM spaxels is their location in the galaxy. Hwang et al. (2019) showed that they occur preferentially at larger radii than the normal spaxels. In Figure 6 we plot all ALM spaxels in the N/O versus O/H plane color-coded by the radius at which the spaxel is located, normalized by the galaxy half-light radius (R/R_e). We see that the ALM spaxels that deviate most strongly from the normal N/O versus O/H relation (e.g., the points lying along the left edge of the distribution in the figure) show a strong trend to be located in the outer regions of their host galaxies (R > 1.5R_e). In the context of the mixing model, these locations in the N/O versus O/H plane would correspond to a more dominant fraction of accreted versus pre-existing gas in these locations. This is qualitatively consistent with a model where gas is accreted from outside the galaxy.

Zoom In Zoom Out Reset image size

The fact that the ALM spaxels have an excess N/O for a given O/H rules out the interpretation that they are simply less chemically evolved regions (in which case they would have lower O/H but still follow the normal N/O versus O/H relation). However, as noted in the introduction, there are alternative ways in which gas in a galaxy could depart from the normal N/O versus O/H relation besides the accretion of metal-poor gas. We briefly discuss the possibilities of these alternatives here.

First, strong winds carrying nitrogen-enriched gas from WR (post-main sequence hot and massive) stars can raise N/O in the ISM (Walsh & Roy 1989; López-Sánchez et al. 2007). We use the high S/N stacked spectra described in the Appendix to search for the two strongest features among the series of weak and broad emission features due to WR stars, at ∼4660 and ∼5810 Å. In Figure 7 we show zoomed-in views of these regions in stellar continuum subtracted stacked spectra of ALM spaxels. In no case do we detect the WR feature of broad emission bumps (for examples of such broad emission bumps, see Figure 20 in Loaiza-Agudelo et al. 2020). Note that Brinchmann et al. (2008) found an enhancement in N/O by ∼0.1 dex in galaxies with detectable WR features compared to those without. This suggests that the ISM can be enriched in nitrogen during a phase in which a strong WR population is present. Since we do not detect a strong WR population in the ALM regions, there is no evidence that WR stars are responsible for the enhanced N/O we see.

Zoom In Zoom Out Reset image size

Second, the release of secondary nitrogen by post-main-sequence intermediate-mass stars (4−7M_⊙, Kobayashi et al. 2020) lags the release of oxygen in more massive stars, since these intermediate-mass stars have lifetimes of ∼10⁸ yr. An excess in N/O could be produced in the case of a strong post-starburst system in which no core-collapse supernovae are releasing newly created oxygen, but the intermediate-mass stars produced 10⁸ yr ago in the burst are finally returning secondary nitrogen (e.g., Schaefer et al. 2020). A selective loss of oxygen via supernovae-driven galactic winds during the burst phase could also contribute to a high N/O value during the post-burst phase (Vincenzo et al. 2016). However, we emphasize that this scenario is inconsistent with the properties of the ALM spaxels, which are shown to be tracing on-going bursts of star formation with a significant population of short-lived massive stars (Hwang et al. 2019).

Finally, an interesting variant of the mixing scenario is one in which metal-rich gas is ejected from the central region of a galaxy as a fountain-flow and mixes with relatively metal-poor gas in the outer disk (Belfiore et al. 2015). While this can produce gas with enhanced N/O, it would also enhance O/H and lead to regions with anomalously high metallicity.

Therefore, we conclude that the accretion of metal-poor gas is responsible for the chemical composition of gas in the ALM regions in our case.

In order to measure the fluxes of weak auroral lines, we must stack spectra of ALM spaxels to obtain an adequate S/N. We select three metallicity bins for stacking defined using previously calculated (12 + log(O/H))_obs from the RS32 calibrator: Z_RS32 < 8.5, 8.5 < Z_RS32 < 8.6, 8.6 < Z_RS32. The spectra to be stacked are taken from unbinned DRP products accessed via marvin (Cherinka et al. 2019). Before stacking, we check both DRP Maskbits and DAP Maskbits from marvin and exclude problematic spaxels with bad quality flags. We then follow Andrews & Martini (2013) to stack our spectra.

We first shift each spectrum to the rest-frame using its redshift and the local line-of-sight emission line velocity obtained from DAP maps. We then interpolate each spectrum in a log-linear way with Δlogλ = 10⁻⁴ Å, between the minimal and maximum wavelength out of all spectra in the corresponding metallicity bin. Finally we normalize the spectra to the mean flux between 4400 and 4450 Å and co-add (average the flux at each wavelength) all the spectra in the metallicity bin.

To analyze the stacked spectra, we use MaNGA DAP software package v2.4.1¹⁸ (Belfiore et al. 2019; Westfall et al. 2019), which automatically fits and subtracts the stellar continuum, providing measured emission line fluxes and corresponding flux errors. We use the MILES-HC template for the stellar continuum fitting and the MASTAR-HC template for emission line fitting, following the same setting in MaNGA MPL-9 data analysis.

The features of interest are the weak auroral lines used to measure T_e: [O iii]4363, [N ii]5755, [O ii]7320, 7331, and the regions between 4600–4700 and 5750–5850 Å where the strongest WR emission features are located. We use the package to fit the stellar continuum and emission lines for only 200 Å excerpts from the complete spectrum, centered on the auroral line(s) of interest. For other stronger emission lines, the results are from fitting the stellar continuum over the entire spectrum. Ideally we would expect S/N to increase by a factor of $\sqrt{N}$ when stacking N individual spectrum together, however, we find in our case that increasing the number of spectra used in stacking could barely increase S/N once the number exceeds ∼1000, i.e., there is irreducible systematic noise in the original spectra that we use for stacking, probably due to imperfect starlight subtraction. Nevertheless, the stacking procedure described above allows us to successfully measure the fluxes of the auroral lines. The extinction is corrected following Section 3.1 above.

In order to determine the electron temperature (T_e), we need to apply corrections based on the electron density. The electron density n_e is estimated from the ratio of [S ii]6716, 6731 lines, using Equation (5) in Yates et al. (2020). We calculate separate electron temperatures for O ii and O iii (T_e(O iii) and T_e(O ii)) by solving Equation (3) in Yates et al. (2020) iteratively. Our calculation process agrees with Yates' statement that convergence is typically achieved within three iterations for T_e(O iii), and five iterations for T_e(O ii). The auroral line [N ii]5755 is too weak to measure even in the stacked spectra, so we assume T_e(N ii) = T_e(O ii). This assumption is supported by the comparisons done by Curti et al. (2017). In estimating the uncertainties of electron temperatures, we assume all uncertainties stem from the flux uncertainties of the auroral lines [O iii]4363, [O ii]7320, 7331 as other stronger emission lines have much smaller flux uncertainties.

The uncertainties in the fluxes of the weak auroral lines are dominated by systematic errors in the subtraction of features in the stellar continuum. Therefore, rather than use the pipeline-estimated errors, we calculate the flux uncertainty in these lines as the standard deviation of the fluxes in the residual (starlight-subtracted) spectrum integrated over wavelength bins with size roughly equal to the FWHM of the corresponding line (10 Å for [O iii]4363, 20 Å for the blended [O ii]7320, 7331). This standard deviation was calculated around each line(s) over the 200 Å fitting window. The uncertainty of the electron temperature is then calculated using a Monte Carlo method where the observed auroral line flux is deviated from its measured value by a Gaussian with σ equal to its flux uncertainty, and an electron temperature is calculated with this new flux. This process is repeated 1000 times and the standard deviation of the resulting distribution is taken as the uncertainty of the measured electron temperature.

Direct measurements of singly and doubly ionized oxygen abundances O⁺/H⁺ and O⁺⁺/H⁺ are calculated using Equations (6) and (7) in Yates et al. (2020), respectively. The uncertainties in oxygen abundances are derived with the same Monte Carlo method described before, propagating the previously found uncertainties in the electron temperatures. The total oxygen abundance is obtained via O/H = O⁺/H⁺ + O⁺⁺/H⁺.

To obtain direct measurements of the N/O ratio, We assume log(N⁺/O⁺) = log(N/O) and calculate it following Equations (3) and (6) in Izotov et al. (2006).

The uncertainties in the direct method metallicities estimated as discussed above is on the order of 0.01 dex. However, we find larger systematic uncertainties in the direct method (12 + log(O/H))_direct. As we already know that increasing the number of raw spectra in the stacking could no longer improve the S/N of the stacked spectrum once the number of raw spectra exceeds ∼1000, we try to make 10 stacked spectra in metallicity bins Z_RS32 < 8.5 (not many spectra fall here), 8.5 < Z_RS32 < 8.6, 8.6 < Z_RS32 < 8.7, each with ∼2000 raw spectra drawn randomly from all spaxels satisfying the criteria described in Section 2.1 in that metallicity bin. We calculate (12 + log(O/H))_direct for the 10 stacked spectra in each of the metallicity bin, and tabulate their mean and standard deviation in Table 2.

We observe large scatter in the measured direct method oxygen abundances in the higher-metallicity bins. While the scatter is significantly smaller in the Z_RS32 < 8.5 bin, 12 + log(O/H) values in this range lie on the flat part of the N/O versus O/H plot and are thus poor diagnostics of the mixing model. Moreover, the standard deviation of oxygen abundances calculated from the stacked spectra is about 0.1 dex in the two higher-metallicity bins. By our definition of ALM spaxels in Section 2.2, the difference in 12 + log(O/H) values between ALM and non-ALM spaxels is only ∼0.111 dex. Given such large systematic uncertainties in the direct method oxygen abundance, we are unable to reliably distinguish ALM and non-ALM spaxels in the N/O versus O/H plane. Although we therefore use the strong line diagnostics in our chemical composition analysis, the high S/N stacked spectra helped us to look for WR star emission features as discussed in Section 4.2.