Present-day mass-metallicity relation for galaxies using a new electron temperature method

Aims. We investigate electron temperature (Te) and gas-phase oxygen abundance (ZTe) measurements for galaxies in the local Universe (z < 0.25). Our sample comprises spectra from a total of 264 emission-line systems, ranging from individual Hii regions to whole galaxies, including 23 composite Hii regions from star-forming main sequence galaxies in the MaNGA survey. Methods. We utilise 130 of these systems with directly measurable Te(Oii) to calibrate a new metallicity-dependent Te(Oiii)–Te(Oii) relation that provides a better representation of our varied dataset than existing relations from the literature. We also provide an alternative Te(Oiii)–Te(Nii) calibration. This new Te method is then used to obtain accurate ZTe estimates and form the mass – metallicity relation (MZR) for a sample of 118 local galaxies. Results. We find that all the Te(Oiii)–Te(Oii) relations considered here systematically under-estimate ZTe for low-ionisation systems by up to 0.6 dex. We determine that this is due to such systems having an intrinsically higher O abundance than O abundance, rendering ZTe estimates based only on [Oiii] lines inaccurate. We therefore provide an empirical correction based on strong emission lines to account for this bias when using our new Te(Oiii)–Te(Oiii) and Te(Oiii)–Te(Nii) relations. This allows for accurate metallicities (1σ = 0.08 dex) to be derived for any low-redshift system with an [Oiii]λ4363 detection, regardless of its physical size or ionisation state. The MZR formed from our dataset is in very good agreement with those formed from direct measurements of metal recombination lines and blue supergiant absorption lines, in contrast to most other Te-based and strong-line-based MZRs. Our new Te method therefore provides an accurate and precise way of obtaining ZTe for a large and diverse range of star-forming systems in the local Universe.


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. 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(Stasińska , 2005Kewley & 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 lowredshift 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 overestimates 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 massmetallicity 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).

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 (z = 0.02). 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 −1 , respectively (Law et al. 2016).
Galaxies were initially selected to have log(M * /h −2 68 M ) < 10.0. 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(Stasińska , 2005Bresolin 2008). This can lead to an under-estimation of the gasphase 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 A107, page 2 of 26 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 1 . 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 −15 erg s −1 cm −2 .
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. 2015Hirschauer et al. , 2018.

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 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 Z Te estimates, whereas those without represent those with semi-direct Z Te estimates. 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.
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 23 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 23 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 A107, page 3 of 26 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.

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) 2 or the NASA-Sloan Atlas v0_1_2 (NSA, Blanton et al. 2011) 3 .
Given that stellar masses are presented in units of h −2 68 M in this work, we have multiplied all SDSS-DR7 masses by the small factor 0.7 2 /0.68 2 = 1.06 and all NSA masses by the larger factor 1 2 /0.68 2 = 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 lowredshift dataset that are present in both these catalogues with a second-order polynomial given by: log(M * /h −2 68 M ) NSA = 2.81707 + 0.542641x + 0.0181975x 2 , (1) in the range 6.0 < x < 10.0, where x = log(M * /h −2 68 M ) SDSS-DR7 . This yields a range of correction factors between 1.01 for highmass 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.

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 R V = 3.1 +2.7 −1.0 , 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 starforming galaxies, k (λ, R V ), and F cor (λ) = F obs (λ) 10 0.176 E(B−V) gas k (λ,R V ) , where R V = 4.05±0.80, 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 3 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.

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. (2013Nicholls et al. ( , 2014a. Those works utilised the MAPPINGS IV photoionisation code  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 1 D 2 and 1 S 0 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 −3 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: 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; where R obs is the reddening-corrected flux ratio of the particular ionic species, and N e is the electron density in cm −3 . The values of the coefficients a, b, c, and d are also dependent on the A107, page 4 of 26 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: 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) 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 −3 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: where g 1 is the statistical weight of the ground state (g 1 = 4 for O + and 9 for O ++ ), α Hβ is the T e (Oiii)-dependent effective emissivity for H β assuming Case B recombination, E 12 = hc/λ avg is the energy difference between the collisionallyexcited 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, Υ 12 is the temperature-dependent net effective collision strength for the transition in question, β = (h 4 /8π 3 m 3 e k) −1/2 = 115885.4 K −1/2 cm −3 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 Υ 12 .
An estimate of the total oxygen abundance, Z Te , can then be obtained by summing the abundances of these two ionic species; 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. 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.

The T e (Oiii)−T e (Oii) relation
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.

A new empirical T e (Oiii)-T e (Oii) 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, A107, page 5 of 26 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 .
Therefore, we propose a functional form for our new Z Tedependent 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 where a is the hyperbolic semi axis. We then fit a to the directlymeasured Z Te values from our dataset, finding the following linear dependence (see Appendix D for details on our fitting methods): 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.
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.
Equations (9) and (10) can be combined to provide the following equivalent expression for Z Te , 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.
A107, page 6 of 26  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 Z Te 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: Z Te −T e (Oiii) and Z Te -T e (Oii) relations for the same dataset. The large spread in the Z Te −T e (Oii) relation leads to quite discrepant fits being returned by our two statistical fitting methods.

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 lowionisation 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". 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.
Physically speaking, systems with low O ++ /O + ostensibly have a larger O + /H + abundance than O ++ /H + abundance,  Fig. 2. Here, each system is coloured by its direct Z Te . Our new empirically-derived T e (Oiii)-T e (Oii) relation (Eq. (9)) is plotted for discrete values of Z Te for comparison. There is a good correspondence seen between the metallicites of the data and the relation across this plane. 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.

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.

[Oii] 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. 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 −3 and O ++ /H + > 4.4 × 10 −4 . For our full dataset, the mean values are 112 cm −3 and 7.0 × 10 −5 , respectively. Also, for their intermediate O ++ /H + estimate of 5.61 × 10 −4 , 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 −3 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)  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 subsample 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.

Dust extinction
Higher metallicity systems are expected to contain more dust, making any differences between the assumed and actual A107, page 8 of 26 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. Kobulnicky et al. (1999) found that T e (Oiii) can be overestimated 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.

Composite Hii regions
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 highresolution 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 subsamples: 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 subsample" 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.

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 emissionline ratio [Oiii] 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.

Correcting the semi-direct Z
where x = log([Oiii]λλ4959, 5007/[Oii]λλ3726, 3729). The corrected semi-direct oxygen abundance is therefore given by, 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.
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 A107, page 9 of 26 A&A 634, A107 (2020) Fig. 8. Comparison between the direct Z Te and the semi-direct Z Te from our new T e (Oiii)-T e (Oii) relation corrected for the O ++ /O + bias using our least-squares-derived f cor (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.  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. 2012or Andrews & Martini 2013 perform relatively well for our composite or integratedgalaxy 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.
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.

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α-fluxweighted 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.
Our new MZR represents an improvement on previous T ebased 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 leastsquares and Bayesian analyses: in the range 5.67 < log(M * ) < 9.87, where M * is in units of h 2 68 M . 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 strongline-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(M * /h 2 68 M ) = 7.51 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.

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.
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 starforming systems (see Sect. 3). We therefore concur with most previous studies that the oxygen abundance in the ISM of lowredshift 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 anticorrelation 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. Figure 13 compares our new MZR (black solid line) with those from other recent studies of T e -based metallicities using different samples.

T e -based MZRs
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 A107, page 11 of 26  14)) and various strong-line-based MZRs presented by Kewley & Ellison (2008, dotted lines). As for many other Z Te studies, we find a lower normalisation for our MZR compared to the strong-linebased relations, although the slope is quite similar among most cases below log(M * /M ) ∼ 10.0. 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 50 < 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 −1 ) < 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 lowredshift 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.

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. (2009Esteban et al. ( , 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 14)), plotted alongside the MZR from oxygen recombination lines (yellow/orange points, Esteban et al. 2009Esteban et al. , 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).
the metallicity at two disc scale lengths from the measured bluesupergiant abundance gradient.
Stellar masses for the Esteban et al. (2009Esteban et al. ( , 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(Esteban et al. , 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. (2009Esteban et al. ( , 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.

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 metallicitydependent 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.  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 oneto-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 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.

Appendix A: An alternative nitrogen-based calibration
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.
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) subsample 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):

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. A107, page 20 of 26

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.
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 . 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. Points with thick black rings represent systems for which we made direct Z Te estimates, while the remaining points have semi-direct Z Te estimates.