© R. M. Yates et al. 2020

1. Introduction

Galaxies are complicated systems. In particular, the bright HII regions in their interstellar medium (ISM) are subject to a number of complex and interrelated astrophysical processes and span a range of sizes, morphologies, luminosities, temperatures, and spatial distributions (e.g. Hodge & Kennicutt 1983; Kennicutt 1988; Osterbrock 1989). This complexity is no less evident when studying the metal content of HII regions. The standard diagnostics used to measure gas-phase metallicities rely on strong, collisionally-excited nebular emission lines which are sensitive not only to metallicity but also a host of other phenomena, such as nebular excitation, shocks, ionsing radiation field strength, gas pressure, electron density, temperature structure, dust content, diffuse ionised gas (DIG) contamination, and the N/O abundance ratio (e.g. Brinchmann et al. 2008; Stasińska 2010; Kewley et al. 2013; Shirazi et al. 2014; Steidel et al. 2014; Krühler et al. 2017; Sanders et al. 2017; Strom et al. 2017; Pilyugin et al. 2018). Even when only considering the local Universe, some strong-line diagnostics are known to be prone to both metallicity-dependent (e.g. Kewley & Dopita 2002; Erb et al. 2006; Yates et al. 2012; Andrews & Martini 2013) and scale-dependent (Krühler et al. 2017) biases. This all means that the true chemical composition of the ISM in nearby galaxies is still not well understood, and warrants continued investigation.

The most direct way to measure the metallicity, or more precisely, the oxygen abundance Z = 12 + log(O/H), in HII regions is via metal recombination lines (RLs). This method allows for an estimate of the abundance of a given ionic species just from measurements of the relevant RLs (e.g. [OII]λ4651 and Hβ, for O⁺⁺) and their effective recombination coefficients, without significant concern for gas temperatures or the properties of the ionising sources (e.g. Osterbrock 1989; Peimbert et al. 1993; Esteban et al. 2002). However, optical metal RLs are extremely weak and require high-resolution spectroscopy of very nearby systems to be detected (e.g. Esteban et al. 2009, 2014). Alternatively, the absorption lines measured in the photospheres of individual blue supergiant stars have been used to obtain similarly direct estimates of gas-phase metallicities (e.g. Urbaneja et al. 2005; Kudritzki et al. 2008). This method relies on the fact that blue supergiants are relatively bright (M_V ∼ −9.5 mag, Bresolin 2003) and relatively young (∼10 Myr old), making their chemical composition both measureable and a fair representation of the recently star-forming gas. This method is, however, also limited to the very nearby Universe currently (z ≲ 0.0025, i.e. ≲10 Mpc). Although, in the era of E-ELT, the use of brighter red supergiant stars (Lardo et al. 2015; Davies et al. 2017) will push the method out to z ∼ 0.025, and super star clusters (Gazak et al. 2014) out to z ∼ 0.2.

A more practical, yet still relatively direct, alternative is the electron-temperature (T_e) method, which relies on measurements of collisionally-excited auroral lines such as [OIII]λ4363 (e.g. Peimbert 1967; Osterbrock 1989; Bresolin et al. 2009). This method works because the T_e of the line-emitting gas is strongly anti-correlated with its metallicity, due to the important role metal ions play in radiative cooling. However, T_e-based metallicities themselves suffer from certain limitations. For example, the auroral lines required are still relatively weak, making measurements difficult to accomplish currently in both more distant (z ≳ 0.3) and more metal-rich (Z ≳ 0.5 Z_⊙) systems (e.g. Bresolin 2008). Additionally, at high metallicities (Z ≳ Z_⊙), saturation of the [OIII]λ4363 line as well as temperature gradients and fluctuations within the HII regions are expected to bias T_e measurements high (and therefore abundance estimates low) (Stasińska 1978, 2005; Kewley & Ellison 2008). Moreover, T_e studies rely on simplified multi-zone models for the temperature structure in ionised nebulae, and often lack information on the temperature of the O⁺ zone, due to the required auroral emission line doublet at [OII]λλ7320,7330 being either too weak for detection or beyond the wavelength coverage of the spectrograph used. This means that O⁺/H⁺, which can be a significant fraction of the total oxygen abundance (see Sect. 5.2), has to be inferred either from a close proxy temperature such as T_e(NII) or from empirical relations linking T_e(OII) to T_e(OIII). This renders such T_e-based metallicity estimates more “semi-direct” than direct.

Therefore, in this work we compile a dataset of 130 low-redshift individual and composite HII regions for which both T_e(OIII) and T_e(OII) are directly measured. These systems are used to investigate the T_e(OIII)–T_e(OII) relation, and derive a new empirical calibration which accounts for the apparent over-estimates in T_e(OII) at low O⁺⁺/O⁺ that we find. This new relation is then used to obtain accurate measurements of Z_Te for a further 134 systems, providing a new insight into the mass – metallicity relation (MZR) of local galaxies.

This paper is organised as follows: in Sect. 2, we outline the new MaNGA sample utilised in this work. In Sect. 3, we assess the possible biases present in our dataset, given its heterogeneous selection. In Sect. 4, we describe how stellar masses, electron temperatures, and oxygen abundances are obtained, as well as how dust corrections are uniformly applied across our dataset. In Sect. 5, we present our analysis of the T_e(OIII)–T_e(OII) relation, including a new metallicity-dependent calibration and an empirical correction for low-O⁺⁺/O⁺ systems. In Sect. 6, we present the MZR, and compare it with those formed from other direct and indirect methods for obtaining metallicity. Finally, in Sect. 7 we provide our conclusions. In Appendix A, we investigate O⁺/H⁺ estimates based on T_e(NII) via the [NII]λ5755 auroral line, which is a possible alternative to the [OII] auroral line quadruplet. In Appendix B, we provide EW(Hα) maps for our MaNGA systems, as well as tables detailing their flux measurements and derived properties. In Appendix C, we provide tables containing the properties we derive from the additional literature samples considered in this work. And in Appendix D, we describe the statistical methods used to fit the key relations presented.

In this work, we make the following distinction between the two ways T_e-based metallicities are obtained: “Direct Z_Te systems” are those for which O⁺⁺/H⁺ and O⁺/H⁺ can be directly determined via measurements of both the [OIII] and [OII] auroral lines (see Sect. 4.3). “Semi-direct Z_Te systems” are those for which only O⁺⁺/H⁺ can be directly determined, and therefore require an assumed relation between T_e(OIII) and T_e(OII) in order to determine O⁺/H⁺ (see Sect. 5).

Throughout, we assume a Chabrier (2003) stellar IMF, and a dimensionless Hubble parameter of h = 0.68 as determined by the Planck Collaboration from combined CMB and lensing data (Planck Collaboration XVI 2014).

2. MaNGA sample

Here, we outline the new MaNGA sample used in this work.

The basis of our dataset is formed of 12 galaxies selected from the Mapping Nearby Galaxies at APO (MaNGA) survey (). MaNGA utilises integral-field units (IFUs) mounted on the Sloan Foundation 2.5 m Telescope at the Apache Point Observatory to obtain fibre-based spatially-resolved spectroscopy of nearby galaxies. The MaNGA sample is taken from the NASA Sloan Atlas catalogue of the SDSS Main Galaxy Legacy Area with a spatial sampling of ∼1−2 kpc and typical resolution of ∼2000 (Bundy et al. 2015).

Our subset of targets was taken from the MaNGA Product Launches-5 (MPL-5), which contains data for a total of 2778 galaxies observed with MaNGA that are fully reduced and vetted as of May 24, 2016 (Law et al. 2016). These data cubes are sky background subtracted using a a super-sampled sky model made from all of the sky fibers, resulting in a typical accuracy of 10% in the skyline regions of the red camera, out to ∼8500 Å, and an accuracy of 10 − 20% at longer wavelengths. This uncertainty was added in quadrature to extracted spectra in regions where skylines were present.

All galaxies are from the main MaNGA survey, as at the time of our data reduction, no systems from the dwarf galaxies ancillary project (Cano Díaz et al., in prep.) had been observed yet. Reduced, calibrated and sky-subtracted data cubes were downloaded from the SDSS Science Archive Server. The median spatial and spectral resolution of these datacubes is 2.54 arcsec (2.60 arcsec for our MaNGA sample) at FWHM and 72 km s⁻¹, respectively (Law et al. 2016).

Galaxies were initially selected to have . This upper mass limit is imposed due to the likely inaccuracy of T_e-based metallicities at super-solar metallicities. In metal-rich HII regions, strong temperature gradients can cause the measured line-ratio temperature of a given species to deviate significantly from the mean ionic temperature to which the oxygen abundance is actually related (Stasińska 1978, 2005; Bresolin 2008). This can lead to an under-estimation of the gas-phase metallicity when using the T_e method by a factor of typically 2 or 3 (Bresolin 2007).

We then carried out two separate analyses on these 12 MaNGA galaxies. The first analysis uses global spectra of each galaxy, only considering spaxels with Hα equivalent widths of EW(Hα) > 30 Å, to minimise contamination from diffuse ionised gas (DIG) emission. The second analysis utilises the IFU capabilities of MaNGA to study distinct regions within these galaxies. We study large Hα-emitting regions which we term “HII blobs” because their effective spatial resolution is in the range 0.8−1.4 kpc (2.3−2.9″), which is several times larger than that of typical HII regions. We select spaxels with Hα equivalent width (EW) of > 50 Å and Hα signal-to-noise ratio (S/N) of > 50, identifying distinct high-EW(Hα) blobs from the resulting maps by eye (see Appendix B), and then extracting spectra from elliptical regions over each blob.

To be able to accurately measure the nebular emission lines without loss of flux from stellar absorption, we removed the stellar continuum from the total emission spectra by fitting each pixel in the MaNGA data cube with the spectral synthesis code STARLIGHT. Subtracting the best-fit stellar component model from the total emission spectrum at each pixel then left us with the gas-only emission spectrum. Emission line fluxes were then measured using Gaussian fits, where in the case of line doublets, the line widths were constrained so that both lines in the doublet had the same velocity width¹. All fits were checked by eye to verify that the procedure was not fitting a Gaussian to noise. In particular, in the case of weak lines, the best-fit line peak (and thus implied redshift) and line widths were compared to the fits to stronger lines. In those cases where the line was not detected or the best-fit Gaussian was fitting noise, the fit was further constrained by fixing the line width and peak position to the best-fit parameters from fits to stronger lines of the same element. The resulting line fluxes are provided in Table B.1 for our HII blob spectra and Table B.2 for our global spectra.

In both analyses, we limit our study to only those systems with S/N([OIII]λ4363) and S/N([OII]λλ7320,7330) ≥ 3. This criterion was applied to ensure that well-constrained T_e-based oxygen abundances could be obtained. For our sample of 23 HII blobs, the mean S/N of the [OIII]λ4363 line is 6.7, and for the [OII]λλ7320,7330 doublet is 23.1. All our MaNGA galaxies exhibit strong Hα lines, with F_obs(Hα)≥1.0 × 10⁻¹⁵ erg s⁻¹ cm⁻².

Our MaNGA sample provides an interesting and relatively new perspective on the low-redshift mass-metallicity relation because all its galaxies lie on or very close to the main sequence of star formation (see Sect. 3), whereas most other studies of individual galaxies with Z_Te ≳ 8.0 are focused on higher-SFR systems (e.g. Izotov et al. 2006; Hirschauer et al. 2015, 2018).

3. Selection effects

Our combined dataset comprises systems of various different physical sizes, from individual HII regions (e.g. the Bresolin et al. 2009 sample), to composite ISM regions (e.g. our MaNGA sample), to integrated galaxy spectra (e.g. the Ly et al. 2016a sample). It is also assembled from various different studies, each with different selection criteria. Therefore, it is important that selection effects are assessed.

One concern is that our T_e measurements could be dependent on the physical size of the system. It has been established that the blending of emission from multiple HII regions of different temperatures within one spectrum can bias estimates of T_e(OII) high when inferring it from T_e(OIII) measurements (Kobulnicky et al. 1999; Pilyugin et al. 2010; Sanders et al. 2017). However, we show in Sect. 5.2.4 that our new T_e(OIII)–T_e(OII) relation performs equally well for both individual and composite HII-region spectra, in part because we have included a range of system sizes in our calibration sample.

Similarly, the possible inclusion of systems with contaminated [OII]λλ7320,7330 auroral lines could affect our results. We have therefore run our analysis on a number of sub-samples, including (a) only systems with high spatial resolution (as described above), (b) only systems with the theoretically expected ratio of the [OII] auroral line doublet, r′ = [OII]λ7320/[OII]λ7330, (c) only systems with similar T_e(OII) and T_e(NII) measurements, (d) using the method from Izotov et al. (2006) for obtaining O⁺/H⁺ that does not require [OII] auroral lines, and (e) discarding T_e(OII) measurements altogether and using T_e(NII) measurements instead. These sub-samples are discussed in Sect. 5.2.2 and Appendix A. We find our results hold for all of these sub-samples, indicating that our new calibration of the T_e(OIII)–T_e(OII) relation (see Sect. 5.1) is not significantly affected by contamination of key auroral lines.

A potential concern for our analysis of the mass-metallicity relation is that datasets such as ours, which require auroral-line detections, can be biased towards starburst galaxies. To assess this, Fig. 1 shows the M_*–SFR relation for our dataset (coloured points) alongside a typical star-forming sample of 109,678 SDSS-DR7 galaxies below z = 0.3 from the study of Yates et al. (2012) (grey contours). To facilitate a fair comparison, only the 55 galaxies from our dataset with a counterpart in the SDSS-DR7 spectroscopic catalogue are considered, and their stellar masses and SFRs are taken from the SDSS-DR7 for this plot. While we do have some systems which are highly star-forming for their mass, there are also a significant number of systems that lie within the 1σ dispersion of the star-forming main sequence of Elbaz et al. (2007) below log(M_*/M_⊙) = 10.0 (∼45% of our SDSS-matched systems). These galaxies are predominantly from our new MaNGA sample, highlighting its important role in this analysis. We therefore conclude that our dataset is relatively representative of the low-redshift star-forming population.

Fig. 1. M*–SFR relation for 55 galaxies in our low-redshift dataset with counterparts in the SDSS-DR7 catalogue. Local “main sequence” relations from Elbaz et al. (2007) and Renzini & Peng (2015) are also plotted for comparison, alongside a star-forming galaxy sample from the SDSS-DR7 (grey contours and points, Yates et al. 2012). Points with black outlines represent systems with direct ZTe estimates, whereas those without represent those with semi-direct ZTe estimates.

Finally, we also checked that we are not preferentially selecting the lowest-metallicity regions within our MaNGA galaxies by requiring S/N([OIII]λ4363) > 3.0. To do this, we measured the metallicity via the R₂₃ strong-line diagnostic for every HII blob within MaNGA systems 8259-9101 and 8459-9102, in order to determine their rough, relative oxygen abundances. Both these systems contain 7 blobs each, allowing a more insightful comparison than for those systems with fewer distinct regions. We use the theoretically-derived R₂₃ diagnostic provided by Kewley & Dopita (2002) and an iterative procedure to obtain both the ionisation parameter, U, and Z_R23 (see Kewley & Dopita 2002, Sect. 6). We find that those HII blobs which make it into our T_e-based sample do not form a special sub-set of the lowest-metallicity regions. In fact, we are able to detect [OIII]λ4363 in some of the most metal-rich regions in these galaxies.

4. Derived properties

4.1. Stellar masses

The majority of the systems in our dataset have stellar masses taken either from the Sloan Digital Sky Survey MPA-JHU data release 7 (SDSS-DR7, Abazajian et al. 2009)² or the NASA-Sloan Atlas v0_1_2 (NSA, Blanton et al. 2011)³.

Given that stellar masses are presented in units of in this work, we have multiplied all SDSS-DR7 masses by the small factor 0.7²/0.68² = 1.06 and all NSA masses by the larger factor 1²/0.68² = 2.16 to maintain consistency throughout.

When stellar masses are available from both catalogues, we adopted those provided by the NSA as they do not suffer from the systematic flux under-estimation for extended systems that SDSS-DR7 stellar masses do, due to improved background subtraction (Blanton et al. 2011). This is particularly important for dwarf galaxies such as those in our study, as these tend to be nearby and therefore more extended on the sky.

The systems with the largest discrepancies in mass estimate between these two catalogues tend to have low specific SFRs (sSFR ≡ SFR/M_*). Galaxies with typical or high SFRs for their mass (such as those in our dataset) tend to have discrepancies of less than 1 order of magnitude. We have therefore corrected all SDSS-DR7 stellar masses used in this work, by fitting the M_*SDSS-DR7 − M_*NSA relation for the 41 galaxies in our low-redshift dataset that are present in both these catalogues with a second-order polynomial given by:

(1)

in the range 6.0 <  x <  10.0, where . This yields a range of correction factors between 1.01 for high-mass systems to 1.12 for low-mass systems.

Masses for the remaining 46% of systems (from NGC 300, SLSN hosts, and the Berg et al. 2012; Ly et al. 2016b samples) are taken directly from the literature.

4.2. Dust corrections

In order to make our dataset as homogeneous as possible, we apply exactly the same reddening corrections to all line fluxes for all systems. For those literature samples where uncorrected line fluxes are not provided, we first re-redden the fluxes by reversing the particular correction used in that work, and then correct all observed fluxes using the following two-step process. Firstly, we correct all lines for Milky Way dust extinction using the Cardelli et al. (1989) extinction law, an extinction factor of , and Milky Way reddening along the line of sight from the Schlafly & Finkbeiner (2011) Galactic reddening map. A fit to the wavelength-dependent uncertainty on the extinction law provided by Cardelli et al. (1989), as well as on their measured value of R_V, is also folded-in to the error propagation for our dust corrections.

Secondly, we correct for internal attenuation within the host system using the Calzetti et al. (2000) attenuation law for star-forming galaxies, , and

(2)

where , and E(B − V)_gas is the colour excess of the ionised gas (see Calzetti 1997). An intrinsic Balmer decrement of 2.86 is assumed, which is suitable for case B recombination in local star-forming galaxies with N_e ∼ 100 cm³ and T_e ∼ 10 000 K (Osterbrock 1989).

We note that this method always returns a non-zero correction to the observed fluxes, even when the internal extinction is determined to be zero, due to the ever-present Milky Way extinction to the redshifted line. Applying uniform dust reddening corrections in this way decreases slightly the scatter in the final electron temperature and oxygen abundance distributions.

4.3. Electron temperatures and oxygen abundances

In order to calculate electron temperatures (T_e) and T_e-based metallicities (Z_Te) for our dataset, we adopt the formalism developed by Nicholls et al. (2013, 2014a,b). Those works utilised the MAPPINGS IV photoionisation code (Dopita et al. 2013) to fit relations between collisionally excited line (CEL) flux ratios from observed spectra to key physical properties of the gas. These relations allow a dependence on the electron density and electron energy distribution, and incorporate updated atomic data, including revised collision strengths for many ionic species. The atomic datasets used in MAPPINGS IV are listed in Table 1 of Nicholls et al. (2013). The collision strengths for the key oxygen ions used in this work are taken from Tayal (2007) for O⁺, and from Aggarwal (1993), Lennon & Burke (1994), Aggarwal & Keenan (1999), and Palay et al. (2012) for O⁺⁺ (with those for the ¹D₂ and ¹S₀ levels from which the λλ4959,5007 and λ4363 lines originate taken from Palay et al. 2012).

The absolute accuracy of collision strength and transition probability estimates is very difficult to constrain, due to the complex calculations and assumptions involved in their determination. Therefore, we have included an additional error on all the electron densities, electron temperatures, and ionic abundances we calculate in this work, to account for the uncertainty in atomic data. This could be particularly important for our T_e(OIII) and O⁺⁺/H⁺ estimates, as the O⁺⁺ collision strengths provided by Palay et al. (2012) are known to return electron temperatures for HII regions that are up to ∼600 K lower than other Breit-Pauli R-matrix based methods from the literature (see Storey et al. 2014; Izotov et al. 2015; Tayal & Zatsarinny 2017).

We account for such atomic data discrepancies by adopting the typical uncertainties found for HII regions with N_e <  1000 cm⁻³ by Juan de Dios & Rodríguez (2017) when comparing 52 different atomic datasets from the literature, including that of Palay et al. (2012). This leads to the following additional errors, which are propagated through all the calculations we make hereafter: σ(N_e) = 0.16 dex, σ[T_e(NII)] = 0.02 dex, σ[T_e(OII)] = 0.02 dex, σ[T_e(OIII)] = 0.025 dex, σ[O⁺/H⁺] = 0.1 dex, σ[O⁺⁺/H⁺] = 0.075 dex. These additional uncertainties increase the typical error on our Z_Te estimates by 0.035 dex.

In this work, we assume thermal equilibrium in the gas (i.e. a Maxwell-Boltzmann electron energy distribution), adopting the following expression for the electron temperature in Kelvin;

(3)

where R_obs is the reddening-corrected flux ratio of the particular ionic species, and N_e is the electron density in cm⁻³. The values of the coefficients a, b, c, and d are also dependent on the ionic species and electron density, and have been calculated by Nicholls et al. (2013, their Tables 4 and 5) by fitting to the output from the MAPPINGS IV code. The flux ratios we consider in this work for R_obs are:

(4)

where the nitrogen ratio, R(NII), is only used for a subset of systems for which the particularly weak [NII]λ5755 auroral line was detected (see Appendix A). We then solve Eq. (3) iteratively to obtain T_e. Convergence is typically achieved within three iterations for T_e(OIII), and five iterations for T_e(OII).

We checked that the electron temperatures obtained from this procedure are very similar to those obtained by numerically solving the complete statistical equilibrium equations using the latest version of the pyneb package (Luridiana et al. 2013), which is a revised version of the nebular/temden routines provided by IRAF (Shaw & Dufour 1995). When assuming the same atomic data, we find the median difference in direct T_e(OIII) is only ∼40 K for our dataset, and the median difference in direct T_e(OII) is only 132 K. This indicates that analytic expressions such as Eq. (3) above (or Eq. (5.4) from Osterbrock & Ferland 2006) are valid in the temperature and density regime studied here.

We did not make any additional corrections to flux ratios to account for Eddington bias (Eddington 1913) in this work, as these are believed to be negligible and depend sensitively on the instrument and type of observation made. However, we note that when applying the corrections to R(OIII) provided by Ly et al. (2016a) from their MMT and Keck data, the mean Z_Te for our sample increases by only 0.026 dex, due to the reduction in the assumed [OIII]λ4363 line strength for systems with S/N([OIII]λ4363) ≲ 7.

An estimate of N_e is obtained from the ratio of [SII] lines, as provided by O’Dell et al. (2013) based on the work of Osterbrock & Ferland (2006);

(5)

Although T_e(OIII) exhibits a negligible dependence on the electron density for the typical range of N_e observed in local HII regions, the MAPPINGS IV code does infer an important effect for T_e(OII), such that electron densities above 100 cm⁻³ would lead to an over-estimate in T_e(OII) by up to 2000 K if not accounted for.

Following Nicholls et al. (2014a), and making the common assumption of a constant collision strength and a uniform electron temperature in each ionic zone, we then obtain singly- and doubly-ionised oxygen abundances as follows:

(6)

(7)

where g₁ is the statistical weight of the ground state (g₁ = 4 for O⁺ and 9 for O⁺⁺), α_Hβ is the T_e(OIII)-dependent effective emissivity for H_β assuming Case B recombination, E₁₂ = hc/λ_avg is the energy difference between the collisionally-excited state and the ground state (where λ_OII = 3727 Å and λ_OIII = 4997 Å are the flux-weighted average wavelengths for the [OII] and [OIII] lines considered here), k is Boltzmann’s constant, Υ₁₂ is the temperature-dependent net effective collision strength for the transition in question, K^−1/2 cm⁻³ s is the constant factor from the collision rate coefficient, h is Planck’s constant, and m_e is the mass of an electron. We refer the reader to Nicholls et al. (2014a, Sect. 4.1) for more details, including the equations used to determine α_Hβ and Υ₁₂.

An estimate of the total oxygen abundance, Z_Te, can then be obtained by summing the abundances of these two ionic species;

(8)

assuming that higher and lower ionised states of oxygen have a negligible contribution (e.g. Stasińska et al. 2012). We note here that this simple addition of the two measured oxygen abundances assumes that either the O⁺ and O⁺⁺ zones are co-spatial, or the amount of H⁺ in each zone is the same.

For many systems, the weak [OII]λλ7320, 7330 line doublet is either not detected or is redward of the wavelength range of the spectrograph used. In such cases, a direct measurement of the [OII] temperature is not possible using Eq. (3), and one must instead rely upon either an alternative ionic temperature (see Appendix A) or an approximation inferred from the measured [OIII] temperature (e.g. Campbell et al. 1986; Garnett 1992; Pagel et al. 1992; Izotov et al. 2006; Pilyugin 2007; Pilyugin et al. 2009; López-Sánchez et al. 2012; Andrews & Martini 2013; Nicholls et al. 2014a). Such T_e(OIII)– T_e(OII) relations are discussed in the following section.

5. The T_e(OIII) − T_e(OII) relation

Figure 2 shows T_e(OIII) versus T_e(OII) for our 130 direct-Z_Te systems. Each system is coloured to indicate the base sample to which it belongs. Six empirically- or theoretically-derived T_e(OIII)–T_e(OII) relations from the literature are also plotted for comparison.

Fig. 2. Te(OIII)–Te(OII) plane containing our 130 systems with direct ZTe estimates [i.e. those for which both Te(OIII) and Te(OII) can be calculated from auroral lines]. Six empirically- or theoretically-derived Te(OIII)–Te(OII) relations from Campbell et al. (1986, black), Pagel et al. (1992, orange), Izotov et al. (2006, blue, for metallicities of 7.2, 7.6, and 8.2), Pilyugin et al. (2009, red), López-Sánchez et al. (2012, grey), and Andrews & Martini (2013, green, as defined in Ly et al. 2016a), are also plotted for comparison. Our dataset reveals a much broader distribution of Te(OII)/Te(OIII) ratios than is expected from the literature relations.

The mean S/N of the auroral lines used here are S/N([OIII]λ4363) = 18.2 and S/N([OII]λλ7320, 7330) = 22.5. The S/N of the nitrogen auroral line for the 53 systems with reported detections is S/N([NII]λ5755) = 8.2 (see Appendix A).

We first note that the distribution of our dataset across the T_e(OIII)–T_e(OII) parameter space is quite broad, including a significant fraction of systems exhibiting lower T_e(OII) than T_e(OIII) (see also Kennicutt et al. 2003). Consequently, there appears to be no clear one-to-one correlation between T_e(OIII) and T_e(OII) in Fig. 2, indicating that none of the common T_e(OIII)–T_e(OII) relations plotted is particularly representative of the true range of T_e(OII)/T_e(OIII) ratios in this dataset. A similar conclusion can be drawn from the samples considered by Izotov et al. (2006, Fig. 4a) and Andrews & Martini (2013, Fig. 3), and is also discussed by Yan (2018) in reference to their Cloudy 17.00 (Ferland et al. 2017) modelling.

Motivated by this issue, in the following sections we develop a new empirical T_e(OIII)–T_e(OII) relation, which allows for a broader range of T_e(OII)/T_e(OIII) ratios.

5.1. A new empirical T_e(O III)–T_e(O II) relation

The form of our new T_e(OIII)–T_e(OII) relation is motivated by two important considerations. First, that electron temperature and oxygen abundance are anti-correlated. Consequently, we would expect systems with T_e(OIII) < T_e(OII) to have O⁺⁺/H⁺ >  O⁺/H⁺, and for their total Z_Te to be relatively insensitive to T_e(OII). This is illustrated by Region A in the schematic in Fig. 3. Likewise, we would expect systems with T_e(OIII) > T_e(OII) to have O⁺⁺/H⁺ <  O⁺/H⁺, and therefore their total Z_Te to be relatively insensitive to T_e(OIII) (see Region B in Fig. 3). Second, that there is an empirical anti-correlation between T_e(OIII) and T_e(OII) at fixed Z_Te for our dataset, (see coloured points in Fig. 5 below). This is also expected theoretically from the equations laid-out in Sect. 4.3, as O⁺/H⁺ ∝ −O⁺⁺/H⁺ at fixed Z_Te. This is in contrast to what is typically assumed for other metallicity-dependent T_e(OIII)–T_e(OII) relations in the literature (e.g. Izotov et al. 2006; Nicholls et al. 2014a). The anti-correlation we find is illustrated by Region C in Fig. 3 by grey lines of constant Z_Te.

Fig. 3. Schematic diagram qualitatively illustrating the expected relation between Te(OIII) and Te(OII) at fixed ZTe in different regions of the Te(OIII)–Te(OII) plane. Grey lines illustrate arbitrary tracks of constant ZTe. Regions A and B denote the approximate regimes where ZTe is likely dominated by only one ionisation state of oxygen (O++/H+ or O+/H+, respectively). Region C denotes the approximate regime where both ionisation states are likely significant, and a uniform anti-correlation between Te(OIII) and Te(OII) at fixed ZTe is expected. We stress that the simplified picture illustrated here may change if, for example, the nebulae are not radiation bounded, or if the density distribution in the O++ and O+ zones are very different.

Therefore, we propose a functional form for our new Z_Te-dependent T_e(OIII)–T_e(OII) relation which broadly embodies these two considerations, and allows for an unrestricted range of possible T_e(OII)/T_e(OIII) ratios. To do this, we adopt the equation of a rectangular hyperbola, centered on 0 Kelvin, given by

(9)

where a is the hyperbolic semi axis. We then fit a to the directly-measured Z_Te values from our dataset, finding the following linear dependence (see Appendix D for details on our fitting methods):

(10)

We find that Z_Te is more tightly correlated with a than with T_e(OIII) or T_e(OII) alone. This is illustrated in Fig. 4, which shows the Z_Te–a relation for our direct dataset in the top panel, and the Z_Te–T_e(OIII) and Z_Te–T_e(OII) relations in the bottom panels. The standard deviation from least-squares fitting is slightly lower for the Z_Te–a relation [with σ(Z_Te) = 0.10, 0.12, and 0.22, respectively], and this relation is marginally favoured over the T_e(OIII) or T_e(OII) relations by our Bayesian analysis, with comparative Bayes factors of 1.2 and 2.3, respectively.

Fig. 4. Top panel: relation between ZTe and the hyperbolic semi-axis a for our direct dataset. Linear fits to the data (Eq. (10)) are shown for our Bayesian analysis (dark grey solid line) and least squares fitting (black solid line). Dashed lines indicate the 1σ scatter around the distribution of ZTe values for the least-squares fit. The grey area around the Bayesian best-fit indicates the range of possible fits considering the full covariance matrix. Bottom panels: ZTe − Te(OIII) and ZTe–Te(OII) relations for the same dataset. The large spread in the ZTe − Te(OII) relation leads to quite discrepant fits being returned by our two statistical fitting methods.

Equation (9) is plotted for discrete values of Z_Te in Fig. 5, alongside our direct Z_Te systems which are coloured by their direct Z_Te. The correspondence between our new relation and the measured Z_Te for our dataset is good across the whole T_e(OIII)–T_e(OII) plane. Equation (9) can then be solved for both Z_Te and T_e(OII) using fixed-point iteration, similar to the approaches used by Izotov et al. (2006), Pilyugin (2007), and Nicholls et al. (2014a), except that we are explicitly accounting for the interdependence of Z_Te and electron temperature here.

Fig. 5. Te(OIII)–Te(OII) plane containing our 130 systems with direct ZTe estimates, as shown in Fig. 2. Here, each system is coloured by its direct ZTe. Our new empirically-derived Te(OIII)–Te(OII) relation (Eq. (9)) is plotted for discrete values of ZTe for comparison. There is a good correspondence seen between the metallicites of the data and the relation across this plane.

Equations (9) and (10) can be combined to provide the following equivalent expression for Z_Te,

(11)

An alternative fit to the Z_Te–a relation, when using T_e(NII) rather than T_e(OII) to determine O⁺, is discussed in Appendix A.

5.2. A semi-direct abundance deficit at low O⁺⁺/O⁺

First, we compare the direct and semi-direct Z_Te estimates for our dataset, using each of the six literature T_e(OIII)–T_e(OII) relations considered here. Significantly, we find that all the literature relations under-estimate Z_Te by up to 0.6 dex for low-ionisation systems. This is illustrated in Fig. 6, where the difference between the semi-direct and direct Z_Te estimates is plotted against O⁺⁺/O⁺. We hereafter refer to this Z_Te under-estimation as the “semi-direct Z_Te deficit”.

Fig. 6. Semi-direct ZTe deficit as a function of O++/O+ when semi-direct ZTe is calculated using each of the six literature Te(OIII)–Te(OII) relations considered in this work. Labels in each panel denote which Te(OIII)–Te(OII) relation is considered. A similar (although often weaker) anti-correlation is also seen for these literature relations.

This deficit at low O⁺⁺/O⁺ is also present for our new T_e(OIII)–T_e(OII) relation, as shown in panel a of Fig. 7. The apparent ubiquity of the semi-direct Z_Te deficit suggests that it is an intrinsic issue for all semi-direct methods that rely on only T_e(OIII) measurements. We have also checked this analysis by using the equations provided by Izotov et al. (2006) to calculate T_e(OIII), O⁺/H⁺, and O⁺⁺/H⁺ (their Eqs. (1)–(5)) rather than Eqs. (3), (6), and (7) above, and find that our results are unchanged.

Fig. 7. Panel a: semi-direct ZTe deficit (as defined in Sect. 5.2) as a function of O++/O+. There is a tight anti-correlation between these two properties at low values of O++/O+. Panel b: semi-direct ZTe deficit as a function of the nebular line ratio [OIII]/[OII]. As in Fig. 4, the Bayesian fit to this relation is given by the solid dark grey line, and the least-squares fit is given by the solid black line. These fits can be used to help correct semi-direct ZTe estimates for the O++/O+ bias seen in panela.

Physically speaking, systems with low O⁺⁺/O⁺ ostensibly have a larger O⁺/H⁺ abundance than O⁺⁺/H⁺ abundance, meaning that singly-ionised oxygen dominates their total oxygen budget. Andrews & Martini (2013) have also determined that a significant fraction of higher-mass galaxies appear to have their overall Z_Te dominated by O⁺/H⁺. They find that log(O⁺⁺/O⁺) typically drops below 0.0 at masses above log(M_*/M_⊙)∼8.2 for stacks of SDSS-DR7 galaxies (their Fig. 5c). Similarly, Curti et al. (2017) found that higher-metallicity galaxies have log(O⁺⁺/O⁺) below 0.0. In what follows, we discuss various possible explanations for the semi-direct Z_Te deficit at low O⁺⁺/O⁺ that we find.

5.2.1. Weak auroral lines

Systematic effects on T_e(OIII) due to inaccurate [OIII]λ4363 measurements are unlikely to play a significant role as we only consider spectra with S/N([OIII]λ4363) ≥ 3.0. Similarly, a systematic over-estimation of [OII]λ7320, 7330 fluxes due to poor line fitting is unlikely, as we have also enforced a S/N lower limit of 3.0 for these lines and find an average S/N([OII]λ7320, 7330) of 22.5. There is also no trend present between S/N([OII]λλ7320, 7330) and the measured semi-direct Z_Te deficit.

5.2.2. [O II] line contamination

The [OII]λλ7320, 7330 auroral line quadruplet lies in a region of the optical spectrum which is populated by a large number of skylines. Additionally, collisional de-excitation, reddening effects, the telluric emission of OH bands, absorption of water bands, and absorption features in the underlying stellar continuum can affect their flux measurements (see e.g. Kennicutt et al. 2003; Pilyugin et al. 2009; Croxall et al. 2015).

When considering our own dataset, we note that collisional de-excitation is already accounted for in the Nicholls et al. (2013) photoionisation models, via an explicit dependence of T_e on electron density. This means that the T_e(OII) and T_e(NII) temperatures we derive for each system are actually quite similar for most of our dataset, as shown in Fig. A.1. We have also been particularly careful to accurately and consistently correct all emission line fluxes for reddening (see Sect. 4.2), and our MaNGA spectra have been corrected for stellar absorption (see Sect. 2). Also, the majority of our dataset lies at redshifts above z ∼ 0.006, meaning the [OII] auroral lines are redshifted away from strong potential contamination from the OH Meinel band emission in Earth’s lower atmosphere.

We also expect contamination from recombination emission to only have a minor effect. This effect in the [OII] auroral lines should be weakest at low O⁺⁺/O⁺, indicating that the semi-direct Z_Te deficit we find at low-O⁺⁺/O⁺ is not caused by such contamination. Liu et al. (2000) detect recombination emission in the planetary nebula NGC 6153. However, the densities and oxygen abundances of NGC 6153 are much higher than for any of the HII regions in our dataset, with N_e >  2000 cm⁻³ and O⁺⁺/H⁺ >  4.4 × 10⁻⁴. For our full dataset, the mean values are 112 cm⁻³ and 7.0 × 10⁻⁵, respectively. Also, for their intermediate O⁺⁺/H⁺ estimate of 5.61 × 10⁻⁴, Liu et al. (2000) state that the contamination from recombination excitation is roughly the same for the [OII]λλ3726, 3729 nebular lines and the [OII]λλ7320, 7330 auroral lines, so that their ratio is relatively unaffected. Similarly, García-Rojas et al. (2018) studied nine planetary nebula with n_e >  3500 cm⁻³ and found T_e(OII) was over-estimated by only a few hundred Kelvin due to recombination excitation. We therefore expect our lower-density HII regions to be even less affected.

To test for any other residual contamination of the [OII]λλ7320, 7330 lines, we have re-run our analysis using the equations provided by Izotov et al. (2006, their Eqs. (3) and (4)), which allow independent estimates of O⁺/H⁺ to be made using either the auroral [OII]λλ7320, 7330 lines or the nebular [OII]λ3727 lines [along with their semi-direct T_e(OIII)–T_e(OII) relation in both cases]. We find that the O⁺/H⁺ estimates obtained from these two methods differ on average by only ∼9% (i.e. ∼0.04 dex), and removing the very few systems with a difference greater than 0.06 dex has no impact on our results. This indicates that the [OII]λλ7320, 7330 lines in our dataset are reliable for obtaining direct estimates of T_e(OII) and O⁺/H⁺.

Alternatively, specific narrow emission (or absorption) features could cause contamination of one of the [OII]λλ7320, 7330 doublets relative to the other, decreasing the sensitivity of our flux estimates. To test if such contamination is significant, we have measured the r′ = [OII]λ7320/[OII]λ7330 ratio for the 102 systems in our dataset where these doublets are resolved. We find a mean r′ of 1.19, which is in good agreement with the expected theoretical value of r′∼1.24 for systems of similar electron density and temperature (Seaton & Osterbrock 1957; De Robertis et al. 1985). However, the scatter in r′ we find is quite large, with σ(r′) = 0.22. We therefore create a sub-sample of 82 systems with r′ values within the typical range accurately measured for nearby HII regions and planetary nebulae: 1.0 <  r′  ≤  1.6 (Kaler et al. 1976; Keenan et al. 1999). We find both our calibration of the T_e(OIII)–T_e(OII) relation and the semi-direct Z_Te deficit at low O⁺⁺/O⁺ are unchanged when using this sub-sample, indicating that neither of these results are driven by inaccurate [OII]λλ7320, 7330 measurements.

Despite these checks all suggesting line contamination does not significantly affect our metallicity calculations, in Appendix A we also provide a separate T_e calibration using nitrogen lines, which allows an estimate of O⁺ to be made without relying on the [OII] auroral lines at all. This alternative [NII]-based calibration has qualitatively the same features as our [OII]-based calibration.

5.2.3. Dust extinction

Higher metallicity systems are expected to contain more dust, making any differences between the assumed and actual attenuation curve more significant when measuring their emission line fluxes. Additionally, the emission lines required to calculate T_e(OII) are widely separated in wavelength, meaning that T_e(OII) would be more sensitive than T_e(OIII) to such dust effects. This could lead to an under-estimation of Z_Te if, for example, galaxy attenuation curves were to systematically flatten with increasing O⁺/H⁺. Such a scenario would imply that it is not the empirical T_e(OIII)–T_e(OII) relations that are inaccurate, but rather the directly measured [OII] temperatures.

However, we find no systematic correlation between the internal extinction, A_V, and the observed semi-direct Z_Te deficit for our systems with direct T_e(OII) measurements. Also, the tight and systematic anti-correlation between the semi-direct Z_Te deficit and O⁺⁺/O⁺ we find suggests that general variations in attenuation curves among systems are not significant here.

5.2.4. Composite H II regions

Kobulnicky et al. (1999) found that T_e(OIII) can be over-estimated by up to ∼1600 K for composite spectra, and consequently that Z_Te can be under-estimated by up to 0.2 dex. Similarly, Pilyugin et al. (2010) demonstrated that composite spectra can return higher T_e(OIII) [and/or lower T_e(OII)] than the mean electron temperature of their constituent HII regions. Though we note that this effect was found to decrease significantly when the composite spectrum contained more than 2 HII regions. Additionally, Sanders et al. (2017) have shown that emission from diffuse ionised gas (DIG) can combine with the above bias to produce a total over-estimation of T_e(OIII) and under-estimation of T_e(OII) of up to ∼2000 K for global galaxy spectra.

When considering this issue for our dataset, we first note that the semi-direct Z_Te deficit we find is present for both composite and individual HII region spectra. For example, the high-resolution spectra analysed by Bresolin et al. (2009) and the CHAOS team should not be prone to the effects attributed to composite spectra, however we find that they also suffer Z_Te discrepancies of up to 0.35 and 0.5 dex, respectively [see panel a of Fig. 7, light blue and yellow points]. Indeed, Pilyugin et al. (2010) also found the same result for some of the high-resolution spectra from the Bresolin et al. (2009) sample (their Fig. 5). Moreover, the composite spectra from our MaNGA samples have been selected to have at least EW(Hα) > 30 Å, which is well above the level expected for DIG regions (e.g. Belfiore et al. 2016), and all the spaxels considered in our global analysis lie within the star-forming or composite regions of the Baldwin et al. (1981, BPT) diagram. Yet some of these systems are also affected by Z_Te discrepancies.

We have further checked the issue of composite-spectra bias by splitting our 130 direct-Z_Te systems into three distinct sub-samples: an “individual HII region sub-sample” containing spectra of 83 very nearby extragalactic HII regions (i.e. Esteban et al. 2009; Bresolin et al. 2009; CHAOS), a “composite HII region sub-sample” containing 27 composite HII regions (i.e. Guseva et al. 2009; MaNGA), and an “integrated galaxy sub-sample” containing integrated spectra from 20 galaxies (i.e. Lee et al. 2004; Izotov et al. 2012; Hirschauer et al. 2015; SLSN hosts). We find that the relations shown in Fig. 7 are very similar for all three of these sub-samples, as are their fits to the Z_Te–a relation. This indicates that our new T_e method is robust to differences in spatial resolution, and that the semi-direct Z_Te deficit we find is a real feature of low-O⁺⁺/O⁺ systems.

5.2.5. Ionisation factor

The final effect we consider here is that of the ionisation state of the gas. This can be defined by the ionisation parameter, U, which is typically approximated by the emission-line ratio [OIII]λλ4959, 5007/[OII]λ3727 or [SIII]λ9067/[SII] λλ6717, 6731 (Kewley & Dopita 2002).

Panel b of Fig. 7 shows that systems with low [OIII]/[OII] indeed tend to have large Z_Te discrepancies when using our new T_e(OIII)–T_e(OII) relation. This relation is not as tight as that for O⁺⁺/O⁺ in panel a, but there is a clear trend of increasing semi-direct Z_Te deficit with decreasing [OIII]/[OII] for all systems below log([OIII]/[OII]) ∼ 0.25. A similar finding was discussed by Kobulnicky et al. (1999), who concluded that semi-direct methods using only [OIII]λ4363 measurements will inevitably under-estimate Z_Te for low-ionisation nebulae, such as those with extended ionised filaments or shells.

We note that inconsistencies in the relationship between [OIII]/[OII] and U could affect the interpretation of Fig. 7 and that a detailed investigation into the link between ionisation state, its proxies using various strong emission lines, and semi-direct Z_Te estimations is beyond the scope of this current work. Therefore, we only conclude here that the estimated Z_Te is particularly sensitive to T_e(OII) for low-O⁺⁺/O⁺ systems because O⁺ dominates their oxygen budget. This physical effect presents a problem for all semi-direct methods that rely only on [OIII]λ4363 measurements, as they do not have a direct handle on the dominant O⁺ ion.

5.3. Correcting the semi-direct Z_Te deficit at low O⁺⁺/O⁺

Given the findings above, we here derive an empirical correction to the semi-direct Z_Te deficit using the easily-observable ratio of oxygen nebular lines, [OIII]λλ4959, 5007/[OII]λλ3726, 3729 (see panel b of Fig. 7). A linear fit to this relation yields a correction factor, f_cor, given by,

(12)

where x = log([OIII]λλ4959, 5007/[OII]λλ3726, 3729). The corrected semi-direct oxygen abundance is therefore given by,

(13)

Figure 8 shows the comparison between direct Z_Te and the corrected semi-direct Z_Te using our new T_e(OIII)–T_e(OII) relation. There is now a much more consistent correspondence between these two methods across a large range of oxygen abundances. However, despite this clear improvement, we do caution that empirical corrections of this nature do not directly relate the physics linking the two variables in question to the discrepancy observed and are, therefore, liable to provide spurious corrections when applied to inappropriate systems.

Fig. 8. Comparison between the direct ZTe and the semi-direct ZTe from our new Te(OIII)–Te(OII) relation corrected for the O++/O+ bias using our least-squares-derived fcor (Eq. (12)). The uncorrected semi-driect estimates for these systems are shown as grey empty circles connected to their corrected values by a grey dashed line.

We also compare the accuracy of our new semi-direct method with those from the literature. We do this by splitting our sample into individual HII region spectra and composite/integrated-galaxy spectra (as defined in Sect. 5.2.4), and plotting the semi-direct Z_Te deficit distributions from each semi-direct method separately for these two sub-samples in Figs. 9 and 10. These figures clearly show that those literature relations calibrated to single HII regions (e.g. Izotov et al. 2006 or Pilyugin et al. 2009) perform relatively well for our individual HII region spectra, but poorly for composite and global spectra. Conversely, those relations calibrated to composite HII regions (e.g. López-Sánchez et al. 2012 or Andrews & Martini 2013) perform relatively well for our composite or integrated-galaxy spectra, but poorly for the highly-resolved systems. Our new relation, on the other hand, calibrated to a varied dataset, works equally well for all systems of any physical size, with a mean semi-direct Z_Te deficit of almost zero in both cases.

Fig. 9. Distribution of ZTe discrepancies for individual HII regions (as defined in Sect. 5.2.4) when using each of the six literature Te(OIII)–Te(OII) relations considered (red histograms). The distribution for our new semi-direct method (dark grey historgram) when corrected for the observed O++/O+ bias (Eq. (12)) is also shown in each panel. The mean of each distribution is shown as a vertical dashed line.

Fig. 10. Same as Fig. 9, but for composite and galaxy-integrated systems, as defined in Sect. 5.2.4.

When using our dataset as a whole, the standard deviation around the peak of the distribution for our new semi-direct method is 0.08 dex, and there are no systems with Z_Te discrepancies greater than −0.35 dex, regardless of their ionisation state.

6. The local MZR

In Fig. 11, we plot the M_*–Z_Te relation (MZR) for the 118 galaxies that make up our dataset. These galaxies contain the 130 individual and composite HII regions discussed in the preceding sections, as well as additional systems for which semi-direct Z_Te estimates can be obtained. For galaxies with electron temperature measurements for more than one HII region, the Hα-flux-weighted mean Z_Te is used, except for our MaNGA galaxies for which the Z_Te obtained from the global spectra are used. Points with thick black outlines in Fig. 11 represent galaxies with direct Z_Te estimates, while the other points represent galaxies with our semi-direct Z_Te estimates.

Fig. 11. Mass-metallicity relation for the 118 galaxies in our dataset. Galaxies are coloured according to the sample from which their fluxes were taken. Points with thick black rings have direct ZTe estimates [i.e. Te(OIII) and Te(OII) independently calculated], while the remaining points have semi-direct ZTe estimates, i.e. using our new semi-direct method, corrected for the observed O++/O+ bias using our least-squares-derived fcor (see Sect. 5.3). The least-squares and Bayesian linear fits to all the points (Eq. (14)) are shown by the solid black line and solid dark grey line, respectively, with the 1σ dispersions given by dashed lines.

Our new MZR represents an improvement on previous T_e-based relations in the literature predominantly because (a) our dataset is not biased to strongly star-forming systems, and (b) we have corrected for the semi-direct Z_Te deficit in low-O⁺⁺/O⁺ systems discussed in Sects. 5.2 and 5.3.

We find the following linear fit to our MZR, using our least-squares and Bayesian analyses:

(14)

in the range 5.67 <  log(M_*) < 9.87, where M_* is in units of . The 1σ dispersion in Z_Te is 0.21 dex for the least-squares fit from residuals. This spread is larger than the 1σ dispersion of ∼0.10 dex obtained from most strong-line-based MZRs (e.g. Tremonti et al. 2004; Yates et al. 2012), suggesting that the local MZR has a wider scatter than usually assumed. The scatter above our relation is partly due to our low-[OIII]/[OII] correction to semi-direct Z_Te estimates (see Sect. 5.3). For example, the Berg et al. (2012) system CGCG035-007A with log and Z_Te = 8.11 has a low log([OIII]/[OII]) ratio of −0.14, which has caused a Z_Te correction upwards by ∼0.31 dex. However, the other systems close to CGCG035-007A on the MZR actually have relatively high log([OIII]/[OII]) ratio, so their high Z_Te estimates are not due to any correction. We find the scatter below the MZR is driven by unusually-high electron temperatures. For example, the low-metallicity outlier, UGC5340A from the Berg et al. (2012) sample (which is believed to be undergoing a merger), has measured T_e(OIII) = 18 830 K and T_e(OII) = 20 359 K, leading to Z_Te = 7.12.

6.1. Comparison with other MZRs

6.1.1. Strong-line MZRs

In Fig. 12, our MZR is shown alongside fits to the MZRs from the various strong-line metallicity diagnostics presented by Kewley & Ellison (2008). Our new T_e-based MZR has a lower normalisation at fixed mass than MZRs based on strong lines, by 0.15 − 0.55 dex at log(M_*/M_⊙)∼9.0, although the slopes are found to be similar.

Fig. 12. Comparison between our new Te-based MZR (solid black line, Eq. (14)) and various strong-line-based MZRs presented by Kewley & Ellison (2008, dotted lines). As for many other ZTe studies, we find a lower normalisation for our MZR compared to the strong-line-based relations, although the slope is quite similar among most cases below log(M*/M⊙)∼10.0.

A lower normalisation for T_e-based MZRs has been seen many times before (e.g. Stasińska 2005; Lee et al. 2006; Andrews & Martini 2013; Ly et al. 2014; Izotov et al. 2015). This is unlikely to be caused by a preferential selection of low-metallicity systems at fixed mass in our case because our dataset contains a significant contribution of typically star-forming systems (see Sect. 3). We therefore concur with most previous studies that the oxygen abundance in the ISM of low-redshift galaxies is lower than indicated by traditional strong-line diagnostics applied to large star-forming galaxy samples such as the SDSS.

Our findings are qualitatively consistent with the fundamental metallicity relation (FMR, Mannucci et al. 2010), although we are unable to draw any firm conclusions on the anti-correlation between Z_Te and SFR at low-mass, due to the relatively low number of systems at fixed mass and metallicity. We do, however, find a clear increase in star-formation rate (SFR) with both Z_Te and M_* for our dataset.

6.1.2. T_e-based MZRs

Figure 13 compares our new MZR (black solid line) with those from other recent studies of T_e-based metallicities using different samples.

Fig. 13. Comparison between our new Te-based MZR Te-based MZR (solid black line, Eq. (14)) and other recent Te-based MZRs from the literature. Fits to the relations from Lee et al. (2006, orange), Andrews & Martini (2013, green), Izotov et al. (2015, light blue), and Ly et al. (2016b, dark blue) are shown. Stellar masses have been corrected for differences in the assumed dimensionless Hubble parameter where possible.

Our MZR appears to be in best agreement with that of Lee et al. (2006), who utilised Spitzer IR photometry to determine stellar masses and a variety of different semi-direct T_e-based metallicities from the literature for their smaller sample of 27 dwarf irregular galaxies. We determine a larger scatter around the MZR for our larger sample of 118 galaxies, and a lower normalisation.

The slope of the MZR derived by Izotov et al. (2015) (light blue) for 3607 star-forming galaxies from the SDSS-DR7 is significantly shallower than the other T_e-based MZRs considered here. This is likely due to the different sample selection criteria applied. Izotov et al. (2015) selected compact (R₅₀ <  6 arcsec) galaxies with high nebular excitation as given by log([OIII]λ5007/Hβ)≳0.5. They also required a detection of the auroral [OIII]λ4363 line from the relatively shallow SDSS spectroscopy in order to estimate T_e(OIII). These criteria all limit the number of higher-metallicity galaxies in their sample at higher mass, effectively selecting systems with similar physical conditions to higher-redshift galaxies. This is likely causing the lower typical metallicities seen in the Izotov et al. (2015) sample above log(M_*/M_⊙)∼8.0, and therefore the shallower slope compared to other works.

The MZR from Andrews & Martini (2013) has a normalisation ∼0.3 dex higher than ours at log(M_*/M_⊙)∼9.0. Their study utilised T_e measurements from stacked SDSS-DR7 spectra, containing around 2000 galaxies in each bin at log(M_*/M_⊙)∼9.0.

Although Fig. 9 shows that the T_e(OIII)–T_e(OII) relation fit to the Andrews & Martini (2013) data returns semi-direct Z_Te estimates around 0.2 dex higher than our relation, we note that this is the case when applied to spectra of individual HII regions. For the more global galaxy spectra that make up the bulk of our MZR, we actually find a relatively close correspondence between the Z_Te estimated via the AM13 relation and our own (see also Fig. 10).

The range of SFRs found in each sample at log(M_*/M_⊙)∼9.0 are also similar, with − 1.0 < log(SFR/M_⊙ yr⁻¹) < 0.0 in both cases. It is therefore unlikely that differences in the T_e(OIII)–T_e(OII) relations or sample selection biases are responsible for the MZR normalisation difference seen here.

Rather, the predominant cause is likely differences in the estimated direct O⁺/H⁺ obtained from each analysis; Andrews & Martini (2013) find a value of 12 + log(O⁺/H⁺)∼8.5 for their log(M_*/M_⊙)∼9.0 stack (their Fig. 5), compared to values between 7.2 and 8.4 for our direct Z_Te systems of a similar mass. O⁺⁺/H⁺ estimates at fixed mass, on the other hand, are more similar between the two studies.

The higher O⁺/H⁺ estimates they obtain could be partly due to the composite spectra issue discussed by Pilyugin et al. (2010) (see Sect. 5.2.4). However this is unlikely to explain the entire difference in Z_Te seen, as stacked spectra should not be more prone to such issues than global spectra of individual galaxies (Pilyugin et al. 2010; Andrews & Martini 2013; Curti et al. 2017). Limitations with the assumption that galaxies of the same stellar mass should have similar metallicities (see Curti et al. 2017), as well as SDSS fibre aperture effects (see Ellison et al. 2008), could also play a role.

Finally, the results presented by Ly et al. (2016b) on the low-redshift MZR from the Metal Abundances Across Cosmic Time (MACT) survey also suggests higher metallicities than most other T_e-based MZRs considered here. Although O⁺⁺/H⁺ abundances were determined in a very similar way to our methodology, the lack of direct O⁺/H⁺ measurements and reliance on the T_e(OIII)–T_e(OII) relation calibrated to the Andrews & Martini (2013) stacked data is likely causing the higher normalisation, as discussed above.

6.1.3. Alternative direct methods

In Fig. 14, the fit to our low-redshift MZR given by Eq. (14) is shown (black line), compared to MZRs formed using the two methods for obtaining gas-phase oxygen abundances that are generally considered to be most accurate. Large yellow points denote 10 nearby (z <  0.023) galaxies for which O/H has been calculated in the brightest HII regions from faint metal recombination lines by Esteban et al. (2009, 2014). Large blue points denote 15 galaxies for which O/H has been calculated from absorption lines in the photospheres of blue supergiant stars by Kudritzki et al. (2016, and references therein), taking the metallicity at two disc scale lengths from the measured blue-supergiant abundance gradient.

Fig. 14. Fit to our Te-based MZR (solid black line, Eq. (14)), plotted alongside the MZR from oxygen recombination lines (yellow/orange points, Esteban et al. 2009, 2014), the MZR from blue supergiant star photospheres (blue points, Kudritzki et al. 2016), and the MZR from the integrated stellar populations of SDSS-DR7 star-forming galaxies (red points, Zahid et al. 2017).

Stellar masses for the Esteban et al. (2009, 2014) samples are taken from the literature (Sick et al. 2014; Lelli et al. 2014; van Dokkum et al. 2014; López-Sánchez 2010; Skibba et al. 2011; Östlin et al. 2003; Woo et al. 2008; and the NSA catalogue). Stellar masses for the Kudritzki et al. (2016) sample are taken from various other works, as listed in Table 4 of Kudritzki et al. (2012). Where possible, we have ensured that these are corrected for any differences in the assumed IMF or cosmology.

A key finding of this study is that there is remarkably good agreement between the MZR formed from our new T_e analysis and those formed from metal recombination lines and blue supergiant absorption lines. There has been long-standing evidence that T_e methods typically under-predict the O⁺⁺/H⁺ abundance by 0.26 − 0.43 dex compared to RL methods for the same individual HII regions (Esteban et al. 2009, 2014). However, it appears from Fig. 14 that, in an overall sense, our T_e-based MZR is consistent with that formed using RLs and stellar absorption-line spectra.

Additionally, Zahid et al. (2017) have also produced an MZR based on stacked absorption-line spectra from the SDSS-DR7, using sequential single-burst (SSB) SPS models. Their luminosity-weighted MZR is also plotted in Fig. 14 as red points. There is also very good agreement between the Zahid et al. (2017) MZR and those of Kudritzki et al. (2016), Esteban et al. (2009, 2014), and this work. This convergence of various different direct metallicity methods on a consistent MZR is a promising sign than the true metal content of nearby galaxies is being correctly probed by our new analysis.

7. Summary and conclusions

Electron temperatures (T_e) and gas-phase oxygen abundances (Z_Te) have been obtained for 264 emission-line systems in the local Universe (z <  0.25). This dataset comprises a mix of individual, composite, and whole-galaxy spectra, belonging to both starbursting galaxies and galaxies on the star-forming main sequence. The 130 systems with direct measurements of both T_e(OIII) and T_e(OII) are utilised to calibrate a new T_e(OIII)–T_e(OII) relation, which can be used to estimate Z_Te using the [OIII]λ4363 auroral line. The resulting mass – metallicity relation (MZR) for 118 low-redshift, star-forming galaxies with 5.5 ≲ log(M_*/M_⊙)≲10.0 is then compared to previous works. Our key findings and conclusions are as follows:

• Due to its hyperbolic functional form, our new metallicity-dependent T_e(OIII)–T_e(OII) relation (Eq. (11)) allows for a wider range of T_e(OII)/T_e(OIII) ratios than previous relations. The semi axis of the hyperbola is tightly constrained by Z_Te (see Sect. 5.1). Both T_e(OII) and Z_Te can therefore be obtained iteratively for any system with a robust [OIII]λ4363 auroral line detection.

• We find that all the literature T_e(OIII)–T_e(OII) relations considered here, as well as our own, underestimate Z_Te for systems with log(O⁺⁺/O⁺) ≲ 0.0 (see Fig. 7a). After investigating the possible causes for this semi-direct Z_Te deficit, we determine that it is most likely due to the physical dominance of O⁺ ions over O⁺⁺ in the HII regions of these systems, making Z_Te estimates difficult to obtain from measurements of T_e(OIII) alone (see Sect. 5.2).

• We therefore provide an empirically-calibrated correction to our semi-direct Z_Te estimates for low-O⁺⁺/O⁺ systems, based on the easily-observable nebular oxygen line ratio [OIII]/[OII]. Our new method can then return accurate Z_Te estimates for systems of either high or low ionisation, regardless of their spatial resolution. Overall, our semi-direct Z_Te estimates are within a standard deviation of 0.08 dex from the directly-measured values, which is comparable to or better than any of the literature T_e(OIII)–T_e(OII) relations considered here (see Figs. 9 and 10).

• The low-redshift MZR formed using our new Z_Te estimates has a similar slope to most strong-line based MZRs but a lower normalisation, as found by most previous studies (see Sect. 6.1.1). The scatter of σ(Z_Te)∼0.21 we find is larger than is typically found when using strong-line diagnostics.

• When comparing to other T_e-based MZRs, we deduce that any differences are mainly due to sample selection biases or differences in the direct determination of O⁺/H⁺, rather than the particular T_e(OIII)–T_e(OII) relation used when obtaining semi-direct Z_Te estimates. The inclusion of many “star-forming main sequence” galaxies in our dataset makes our MZR more representative of the typically-star-forming population at low redshift.

• Encouragingly, our new T_e-based MZR is in very good agreement with the MZRs obtained via direct metallicity measurements from metal recombination lines or blue supergiant absorption lines (see Sect. 6.1.3). This is a strong indication that our study is accurately probing the true range of metallicities present in the star-forming galaxy population at low redshift.

In a follow-up paper, we will compare our new T_e-based MZR with MZRs derived from alternative direct methods at higher redshifts, in order to ascertain the true evolution of the gas-phase metallicity in galaxies over cosmic time.

Acknowledgments

The authors would like to thank the anonymous referee for very helpful comments and suggestions, as well as Danielle Berg, Fabio Bresolin, I-Ting Ho, Rolf-Peter Kudritzki, Guinevere Kauffmann, Thomas Krühler, Brent Miszalski, Gwen Rudie, Alice Shapley, Martin Yates, and Jabran Zahid for valuable discussions during the undertaking of this work. We would also like to thank Ricardo Amorín, Fabio Bresolin, Alec Hirschauer, Janice Lee, Matt Nicholl, and John Salzer for providing additional data and guidance, and Christophe Morisset for help with running the pyneb package. This research was partly supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. The authors would also like to acknowledge the TOPCAT interactive graphical viewer and editor (Taylor 2005) which was used for quick analysis and visualisation of our tabulated data. R. M. Y., T.-W. C., and P. W. acknowledge the support through the Sofia Kovalevskaja Award to P. Schady from the Alexander von Humboldt Foundation of Germany.

References

Appendix A: An alternative nitrogen-based calibration

In addition to the calibration of our new T_e(OIII)–T_e(OII) relation presented in Sect. 5.1, we also present here a complimentary calibration to the T_e(OIII)–T_e(NII) relation, using a subset of 53 systems (hereafter, the “nitrogen sub-sample”) for which measurements of both the [OIII]λ4363 and [NII]λ5755 auroral lines are available.

[NII]λ5755 is a particularly weak auroral line, with a mean S/N of only 8.2 for our nitrogen sub-sample, compared to S/N([OIII]λ4363) = 17.4 and S/N([OII]λλ7320, 7330) = 27.3 for the same systems. Nonetheless, due to the similar ionisation energies required for [OII] and [NII], both are expected to trace the same ionisation zone within an HII region, and return similar electron temperatures. Therefore, T_e(NII) can be a useful alternative to T_e(OII) when estimating O⁺, in cases where [OII]λλ7320, 7330 is not available or believed to be contaminated (see Sect. 5.2.2).

Figure A.1 shows the R(NII)–R(OII) relation and T_e(NII)–T_e(OII) relation for our nitrogen sub-sample. With the exception of a few systems with unexpectedly high R(OII) and T_e(OII) (outlined in grey), we find a good linear relationship between these line ratios and temperatures, with a median difference between T_e(OII) and T_e(NII) of only 165 K. In particular, the good one-to-one correspondence between T_e(OII) and T_e(NII) for the CHAOS sample (blue points) is a clear improvement on that reported by the original works (Berg et al. 2015; Croxall et al. 2016). This is because we estimate higher T_e(NII) at fixed T_e(OII) than those studies, bringing these two temperatures more in line with each other.

Fig. A.1. Top panel: relation between R(NII) = [NII]λ5755/(1.33⋅[NII]λ6584) and R(OII) = [OII]λλ7320, 7330/[OII]λλ3726, 3729, for our nitrogen sub-sample of 53 systems. Bottom panel: relation between Te(NII) and Te(OII) for the same nitrogen sub-sample. These two temperatures are expected to be roughly the same for a given HII region, as indicated by the line of equality (dashed line). Points with grey outlines (in both panels) are outliers with T(OII) > T(NII) + 2000 K, or equivalently, R(OII) > R(NII) + 0.0175.

The T_e(OIII)–T_e(NII) relation is shown in the top panel of Fig. A.2. There are fewer systems with large differences between these temperatures than seen for the T_e(OIII)–T_e(OII) relation of our full calibration sample (see Fig. 2). The lack of systems in the upper part of this plane (i.e. in region A of the schematic shown in Fig. 3) is mainly due to measured T_e(NII) being lower than T_e(OII) for those systems (outlined in grey in Fig. A.1). The lack of systems in the lower part of this plane (i.e. region B in Fig. 3) is chiefly a selection bias: systems with higher metallicity have particularly weak [NII]λ5755, and are therefore typically removed from samples selected on this line.

Fig. A.2. Top panel: Te(OIII)–Te(NII) relation for our nitrogen sub-sample of 53 systems (points). Fits to the Te(OIII)–Te(OII) relation from the literature are also shown. Centre panel: relation between direct ZTe and the hyperbolic semi-axis, a. Red lines indicate the least-squares fit to our full calibration sample, using Te(OII) (see Eq. (10)). Black lines indicate the least-squares fit to our nitrogen sub-sample when using Te(NII). Bottom panel: relation between the semi-direct ZTe deficit and [OIII]/[OII]. Again, red lines denote our least-squares [OII]-based fit (see Fig. 7), black lines denote our least-squares alternative [NII]-based fit.

The centre panel of Fig. A.2 shows the Z_Te–a relation. The red lines represent the least-squares fit for our full calibration sample, as shown in Fig. 4. We find there is very little difference in the fit obtained when removing those systems that have T(OII) > T(NII) + 2000 K (shown with grey outlines in Fig. A.1). This further demonstrates that our new calibration performs well for systems on and off the expected relations for individual HII regions.

The black lines in the centre panel of Fig. A.2 represent the least-squares fit to our nitrogen sub-sample when using T_e(NII) to obtain O⁺/H⁺. This relation is steeper than that using T_e(OII) to obtain O⁺/H⁺. The steepness is driven by the small number of systems with a ≳ 17 000, which have higher measured T_e(NII) than T_e(OII). It is important to note that all these systems have log(O⁺⁺/O⁺) > 0.5, so are relatively insensitive to how O⁺ is estimated anyway. The majority of the T_e(NII) sub-sample is still well fit by our [OII]-based calibration. Also, the semi-direct Z_Te deficit correction we derive is very similar for both the full calibration sample and nitrogen sub-sample (bottom panel). Nonetheless, we still provide the revised least-squares fits to our hyperbolic semi axis, a, and correction factor, f_cor, when using T_e(NII):

(A.1)

(A.2)

As with our [OII]-based calibration, our [NII]-based calibration returns relatively accurate Z_Te estimates for the full range of systems available (as discussed in Sect. 5.3).

Appendix B: MaNGA EW(Hα) maps and tables

In Fig. B.1, we present the Hα EW maps for our sample of MaNGA galaxies (see Sect. 2). White ellipses signify HII blobs with S/N[OIII]λ4363 ≥ 3.0. All spaxels within an ellipse with EW(Hα) > 50 Å and S/N(Hα) > 50 are used to determine emission line fluxes. Table B.1 provides the emission line intensities measured for each of the MaNGA HII blobs in our sample, along with their derived electron temperatures and oxygen abundances. Table B.2 provides similar data from our global MaNGA galaxy spectra, and includes information on their position, stellar mass, extinction, and electron density.

Fig. B.1. EW(Hα) maps for the 12 galaxies from our MaNGA sample containing HII blobs with auroral line detections (see Sect. 2). HII blobs are shown and numbered in each panel. White ellipses signify HII blobs with S/N[OIII]λ4363 ≥ 3.0, which are used for our Te analysis in Sect. 4.3.

Appendix C: Literature emission-line samples

Here, we provide tables listing the electron temperatures and oxygen abundances measured for direct systems (Table C.1) and semi-direct systems (Table C.2) from the literature samples utilised in this work.

We also systematically compare the Z_Te estimates we obtain following the procedures described in Sects. 4.3 and 5.1 with those provided by the original papers. This comparison is illustrated in Fig. C.1. We find a median difference in Z_Te of only 0.026 dex on a system-to-system basis.

Fig. C.1. Comparison of the Te-based oxygen abundances obtained from our study and those of previous studies for the same systems. Points with thick black rings represent systems for which we made direct ZTe estimates, while the remaining points have semi-direct ZTe estimates.

The values of Z_Te we obtain for some of the Hirschauer et al. (2015) sample when using the [OII] fluxes provided are significantly higher than those obtained when assuming only [OIII] lines and any T_e(OIII)–T_e(OII) relation, with a mean increase of +0.27 dex. This is due to the apparently dominant contribution of O⁺/H⁺ to the total oxygen abundance in these systems (see Sect. 5.2), which we consider to be robust given the relatively high S/N of the [OII] lines, at least compared to that of [OIII]λ4363.

The few metal-rich HII regions from the CHAOS sample with discrepancies greater than ∼0.2 dex are exclusively from NGC 628 (Berg et al. 2015). These discrepancies could be due to the use of T_e(NII) to estimate O⁺/H⁺ in the original work, which they find is systematically lower than their measured T_e(OII), with little correlation between the two temperatures. This issue, and a similar issue between T_e(OIII) and T_e(SIII) also reported by Berg et al. (2015), is not seen in the later CHAOS data for NGC 5457 (Croxall et al. 2016), whose Z_Te estimates agree more closely with ours.

Appendix D: Linear fitting methods

We adopt two distinct methods when making linear fits to the three key relations discussed in this work, namely, the Z_Te–a relation (Eq. (10)), our [OIII]/[OII]-based semi-direct Z_Te correction (Eq. (12)), and the galaxy MZR (Eq. (14)). An overview of these two fitting methods is provided below.

D.1. Least-squares fitting

The first method we adopt is the linear least-squares approximation provided by the IDL routine MPFITEXY (Williams et al. 2010), which depends on the MPFIT package (Markwardt 2009). This routine utilises the χ² minimisation technique to obtain a straight-line best fit to the data, accounting for errors in both the ordinate and abscissa axes and a determination of the intrinsic scatter which is adjusted to ensure the reduced χ² ∼ 1 using the method outlined in Bedregal et al. (2006).

D.2. Nested sampling fitting

In addition to the linear least-squares fitting method outlined above, we also applied a simple Bayesian fit using nested sampling. For this we used the dynesty package (Speagle 2019). For deriving the posterior we used dynamic nested sampling (Higson et al. 2019).

D.2.1 Linear models

For the comparison between the Z_Te–a relation and the Z_Te–T_e(OIII) and Z_Te–T_e(OII) relations (Sect. 5.1), we use the same prior in all cases. We note that the derived evidence depends very strongly on the prior, and that improper priors can lead to inaccurate results.

In general, each relation combines one or more observable X with the desired quantity Z, such that for each observable we have a dataset of X_i and Z_i.

We assume that each measurement, X_i, obs, is drawn from a normal distribution centered at the value of the model function, f_x, for this measurement pair with an effective standard deviation value. This gives,

(D.1)

The effective standard deviation σ_{Z, i, eff} is defined as follows and includes the standard deviation of both observables as well as an additional scattering term f,

(D.2)

where X_i, true is the true value of X_i, which X_i, obs is sampled from (conceptually). Since sampling an N + 3 dimensional cube is quite expensive, we approximated the partial derivative using the observed value. Since the range for values for Z are typically quite small, other errors (systematic uncertainties, for instance) will dominate over the approximation error. In the case where ∂f_x/∂X_i, true does not depend on X_i, true, the solution is exact (as is the case for the simple linear relations).

For most cases in the paper we use the following model function:

(D.3)

This is used for the various temperature – metallicity fits, as well as the MZR fit. We use weak uninformative priors for the parameters:

(D.4)

(D.5)

(D.6)

When using the information from the posterior, one has to keep in mind the whole covariance matrix,

(D.7)

While the off-axis elements between most parameters are mostly zero, the slope α and intercept β are highly correlated.

D.2.2 Z_Te correction factor

To fit the semi-direct Z_Te correction factor discussed in Sect. 5.3, we used a fairly simple linear model that assumes all measurements below a certain critical [OIII]/[OII] ratio, x_c∼-pagination 𝒰(−0.5, 1.0) [dex], under-estimate the metallicity by a [OIII]/[OII]-dependent factor. The location of the cutoff x_c was left as a free parameter with a non-informative prior.

We assumed the following relation for this correction factor:

(D.8)

which is equivalent to Eq. (12).

All Tables

All Figures

Fig. 1. M*–SFR relation for 55 galaxies in our low-redshift dataset with counterparts in the SDSS-DR7 catalogue. Local “main sequence” relations from Elbaz et al. (2007) and Renzini & Peng (2015) are also plotted for comparison, alongside a star-forming galaxy sample from the SDSS-DR7 (grey contours and points, Yates et al. 2012). Points with black outlines represent systems with direct ZTe estimates, whereas those without represent those with semi-direct ZTe estimates.

In the text

Fig. 2. Te(OIII)–Te(OII) plane containing our 130 systems with direct ZTe estimates [i.e. those for which both Te(OIII) and Te(OII) can be calculated from auroral lines]. Six empirically- or theoretically-derived Te(OIII)–Te(OII) relations from Campbell et al. (1986, black), Pagel et al. (1992, orange), Izotov et al. (2006, blue, for metallicities of 7.2, 7.6, and 8.2), Pilyugin et al. (2009, red), López-Sánchez et al. (2012, grey), and Andrews & Martini (2013, green, as defined in Ly et al. 2016a), are also plotted for comparison. Our dataset reveals a much broader distribution of Te(OII)/Te(OIII) ratios than is expected from the literature relations.

In the text

Fig. 3. Schematic diagram qualitatively illustrating the expected relation between Te(OIII) and Te(OII) at fixed ZTe in different regions of the Te(OIII)–Te(OII) plane. Grey lines illustrate arbitrary tracks of constant ZTe. Regions A and B denote the approximate regimes where ZTe is likely dominated by only one ionisation state of oxygen (O++/H+ or O+/H+, respectively). Region C denotes the approximate regime where both ionisation states are likely significant, and a uniform anti-correlation between Te(OIII) and Te(OII) at fixed ZTe is expected. We stress that the simplified picture illustrated here may change if, for example, the nebulae are not radiation bounded, or if the density distribution in the O++ and O+ zones are very different.

In the text

Fig. 4. Top panel: relation between ZTe and the hyperbolic semi-axis a for our direct dataset. Linear fits to the data (Eq. (10)) are shown for our Bayesian analysis (dark grey solid line) and least squares fitting (black solid line). Dashed lines indicate the 1σ scatter around the distribution of ZTe values for the least-squares fit. The grey area around the Bayesian best-fit indicates the range of possible fits considering the full covariance matrix. Bottom panels: ZTe − Te(OIII) and ZTe–Te(OII) relations for the same dataset. The large spread in the ZTe − Te(OII) relation leads to quite discrepant fits being returned by our two statistical fitting methods.

In the text

	Fig. 5. Te(OIII)–Te(OII) plane containing our 130 systems with direct ZTe estimates, as shown in Fig. 2. Here, each system is coloured by its direct ZTe. Our new empirically-derived Te(OIII)–Te(OII) relation (Eq. (9)) is plotted for discrete values of ZTe for comparison. There is a good correspondence seen between the metallicites of the data and the relation across this plane.
In the text

	Fig. 6. Semi-direct ZTe deficit as a function of O++/O+ when semi-direct ZTe is calculated using each of the six literature Te(OIII)–Te(OII) relations considered in this work. Labels in each panel denote which Te(OIII)–Te(OII) relation is considered. A similar (although often weaker) anti-correlation is also seen for these literature relations.
In the text

Fig. 7. Panel a: semi-direct ZTe deficit (as defined in Sect. 5.2) as a function of O++/O+. There is a tight anti-correlation between these two properties at low values of O++/O+. Panel b: semi-direct ZTe deficit as a function of the nebular line ratio [OIII]/[OII]. As in Fig. 4, the Bayesian fit to this relation is given by the solid dark grey line, and the least-squares fit is given by the solid black line. These fits can be used to help correct semi-direct ZTe estimates for the O++/O+ bias seen in panela.

In the text

	Fig. 8. Comparison between the direct ZTe and the semi-direct ZTe from our new Te(OIII)–Te(OII) relation corrected for the O++/O+ bias using our least-squares-derived fcor (Eq. (12)). The uncorrected semi-driect estimates for these systems are shown as grey empty circles connected to their corrected values by a grey dashed line.
In the text

Fig. 9. Distribution of ZTe discrepancies for individual HII regions (as defined in Sect. 5.2.4) when using each of the six literature Te(OIII)–Te(OII) relations considered (red histograms). The distribution for our new semi-direct method (dark grey historgram) when corrected for the observed O++/O+ bias (Eq. (12)) is also shown in each panel. The mean of each distribution is shown as a vertical dashed line.

In the text

	Fig. 10. Same as Fig. 9, but for composite and galaxy-integrated systems, as defined in Sect. 5.2.4.
In the text

Fig. 11. Mass-metallicity relation for the 118 galaxies in our dataset. Galaxies are coloured according to the sample from which their fluxes were taken. Points with thick black rings have direct ZTe estimates [i.e. Te(OIII) and Te(OII) independently calculated], while the remaining points have semi-direct ZTe estimates, i.e. using our new semi-direct method, corrected for the observed O++/O+ bias using our least-squares-derived fcor (see Sect. 5.3). The least-squares and Bayesian linear fits to all the points (Eq. (14)) are shown by the solid black line and solid dark grey line, respectively, with the 1σ dispersions given by dashed lines.

In the text

	Fig. 12. Comparison between our new Te-based MZR (solid black line, Eq. (14)) and various strong-line-based MZRs presented by Kewley & Ellison (2008, dotted lines). As for many other ZTe studies, we find a lower normalisation for our MZR compared to the strong-line-based relations, although the slope is quite similar among most cases below log(M*/M⊙)∼10.0.
In the text

Fig. 13. Comparison between our new Te-based MZR Te-based MZR (solid black line, Eq. (14)) and other recent Te-based MZRs from the literature. Fits to the relations from Lee et al. (2006, orange), Andrews & Martini (2013, green), Izotov et al. (2015, light blue), and Ly et al. (2016b, dark blue) are shown. Stellar masses have been corrected for differences in the assumed dimensionless Hubble parameter where possible.

In the text

	Fig. 14. Fit to our Te-based MZR (solid black line, Eq. (14)), plotted alongside the MZR from oxygen recombination lines (yellow/orange points, Esteban et al. 2009, 2014), the MZR from blue supergiant star photospheres (blue points, Kudritzki et al. 2016), and the MZR from the integrated stellar populations of SDSS-DR7 star-forming galaxies (red points, Zahid et al. 2017).
In the text

Fig. A.1. Top panel: relation between R(NII) = [NII]λ5755/(1.33⋅[NII]λ6584) and R(OII) = [OII]λλ7320, 7330/[OII]λλ3726, 3729, for our nitrogen sub-sample of 53 systems. Bottom panel: relation between Te(NII) and Te(OII) for the same nitrogen sub-sample. These two temperatures are expected to be roughly the same for a given HII region, as indicated by the line of equality (dashed line). Points with grey outlines (in both panels) are outliers with T(OII) > T(NII) + 2000 K, or equivalently, R(OII) > R(NII) + 0.0175.

In the text

Fig. A.2. Top panel: Te(OIII)–Te(NII) relation for our nitrogen sub-sample of 53 systems (points). Fits to the Te(OIII)–Te(OII) relation from the literature are also shown. Centre panel: relation between direct ZTe and the hyperbolic semi-axis, a. Red lines indicate the least-squares fit to our full calibration sample, using Te(OII) (see Eq. (10)). Black lines indicate the least-squares fit to our nitrogen sub-sample when using Te(NII). Bottom panel: relation between the semi-direct ZTe deficit and [OIII]/[OII]. Again, red lines denote our least-squares [OII]-based fit (see Fig. 7), black lines denote our least-squares alternative [NII]-based fit.

In the text

	Fig. B.1. EW(Hα) maps for the 12 galaxies from our MaNGA sample containing HII blobs with auroral line detections (see Sect. 2). HII blobs are shown and numbered in each panel. White ellipses signify HII blobs with S/N[OIII]λ4363 ≥ 3.0, which are used for our Te analysis in Sect. 4.3.
In the text

	Fig. C.1. Comparison of the Te-based oxygen abundances obtained from our study and those of previous studies for the same systems. Points with thick black rings represent systems for which we made direct ZTe estimates, while the remaining points have semi-direct ZTe estimates.
In the text