THE METAL ABUNDANCES ACROSS COSMIC TIME (${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$) SURVEY. II. EVOLUTION OF THE MASS–METALLICITY RELATION OVER 8 BILLION YEARS, USING [O iii] λ4363 Å BASED METALLICITIES

To address the lack of [O iii] λ4363 measurements at higher redshifts, and outstanding issues with gas metallicity calibrations for higher redshift galaxies, we conducted a spectroscopic survey called "Metal Abundances across Cosmic Time" (${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$; Ly et al. 2016, hereafter Paper I) to obtain deep (2–12 hr) rest-frame optical spectra of $z\lesssim 1$ star-forming galaxies with Keck and MMT. The primary goal of the survey was to obtain reliable measurements of the gas-phase metallicity and other physical properties of the interstellar medium (ISM) in galaxies, such as the SFR, gas density, ionization parameter, dust content, and the source of photoionizing radiation (star formation and/or active galactic nucleus, AGN). ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ is unique among previous spectroscopic surveys because it is the first to use the T_e method to measure the evolution of the M_⋆–Z relation over $\approx 8$ billion years. In addition, the galaxy sample of ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ encompasses nearly 3 dex in stellar mass, including dwarfs as low as ${M}_{\star }\sim 3\times {10}^{6}$ M_☉ and 3 × 10⁷ M_☉ at $z\sim 0.1$ and $z\sim 1$, respectively. The ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ survey targeted ≈1900 galaxies in the Subaru Deep Field (SDF; Kashikawa et al. 2004) that have excess flux in narrowband and/or intermediate-band filters, which is now understood to be produced by nebular emission lines from star formation or AGNs (e.g., Ly et al. 2007, 2011, and references therein).

In this paper, Paper II, we focus on the first results from 66 galaxies with at least S/N = 3⁵ detections of [O iii] λ4363 at z = 0.05–0.95 (average of z = 0.53 ± 0.25; median of 0.48), and robust [O iii] λ4363 upper limits for 98 galaxies at z = 0.04–0.96 (average of $z=0.52\pm 0.23;$ median of 0.48). We refer to the collective of these galaxies as the "[O iii] λ4363-detected and [O iii] λ4363-non-detected samples." For the [O iii] λ4363-non-detected galaxies, we require an [O iii]λ5007 detection that is at S/N ≳ 100 and S/N < 3 for [O iii] λ4363. This work expands on our previous sample of spectroscopic detections of [O iii] λ4363 (Ly et al. 2014) by more than threefold. In a forthcoming paper, we will use our sample with [O iii] λ4363 measurements to recalibrate the strong-line metallicity diagnostics for these galaxies at $z\approx 0.5$.

We refer readers to Paper I for more details on the ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ survey and our primary sample for Paper II. Specifically, Section 2 in Paper I describes the full galaxy sample and optical spectroscopy, Section 3 in Paper I describes the [O iii] λ4363-detected and [O iii] λ4363-non-detected sample selection, and Section 4 in Paper I describes the interstellar (i.e., T_e-based metallicity, dust attenuation) and stellar properties (i.e., SFR, stellar mass) of [O iii] λ4363-detected and [O iii] λ4363-non-detected galaxies. The outline of Paper II is as follows.

In Section 2, we discuss the identification of a small number of AGNs or low-ionization nuclear emitting regions (LINERs; Heckman 1980) that contaminate our galaxy sample. In Section 3, we present our five main results: (1) a large sample of extremely metal-poor galaxies at $z\gtrsim 0.1$, (2) comparison of our samples against other star-forming galaxies on the M_⋆–SFR projection, (3) the similarity of these metal-poor galaxies to typical star-forming galaxies at high-z, (4) the evolution of the T_e-based M_⋆–Z relation, and (5) the secondary dependence of the M_⋆–Z relation on SFR. In Section 4, we compare our M_⋆–Z relation against predictions from theoretical and numerical simulations, discuss the selection function of our survey, and compare our survey to previous T_e-based studies. We summarize results in Section 5.

Throughout this paper, we adopt a flat cosmology with ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$, ${{\rm{\Omega }}}_{M}=0.3$, and ${H}_{0}=70$ km s⁻¹ Mpc⁻¹. Magnitudes are reported on the AB system (Oke 1974). For reference, we adopt 12 + $\mathrm{log}({\rm{O/H}})$ ${}_{\odot }$ = 8.69 (Allende Prieto et al. 2001) as solar metallicity, Z_☉. Unless otherwise indicated, we report 68% confidence measurement uncertainties, and "[O iii]" alone refers to the 5007 Å emission line.

We identified a total of 16 extremely metal-poor galaxies with 12 + $\mathrm{log}({\rm{O/H}})$ $\,\leqslant$ 7.69 (i.e., less than 10% of solar). This is the largest extremely metal-poor galaxy sample at $z\gtrsim 0.1$. Keck06 is our most metal-poor galaxy with $12\,+\,\mathrm{log}({\rm{O/H}})={7.23}_{-0.14}^{+0.11}$ (3% of solar metallicity). This is similar to I Zw 18, which is the most metal-deficient galaxy known in the local universe.

We find that 24% of our [O iii] λ4363-detected galaxies are extremely metal-poor; this is far higher than the 4% of [O iii] λ4363-detected galaxies in SDSS that are extremely metal-poor (Izotov et al. 2006). This is presumably attributable to a combination of redshift evolution (lower metallicity toward higher redshift; see Section 3.4) and selection effects, because our sample is focused on lower-mass galaxies (≲10⁹ M_☉) that tend to have lower metallicity. If the extremely metal-poor galaxy fraction increases toward even lower masses, it is possible that a substantial minority of local galaxies—by number—are extremely metal-poor, even though their total mass is only a small fraction of the current total stellar mass in the universe. We suggest that future selections of extremely metal-poor galaxies should either use narrowband imaging or grism spectroscopy. This is more efficient observationally than a brute-force approach within a magnitude-limited survey. For example, the DEEP2 survey (Newman et al. 2013), which targeted ${R}_{{\rm{AB}}}\lesssim 24$ galaxies, has identified only two extremely metal-poor galaxies at $z\sim 0.8$ from a sample of 28 [O iii] λ4363-detected galaxies (Ly et al. 2015).

In Figure 2, we compare our dust-corrected instantaneous SFRs from Hα or Hβ luminosities against stellar masses determined from spectral energy distribution (SED) fitting, to locate our galaxies on the M_⋆–SFR relation and to compare against other star-forming galaxies at $z\lesssim 1$. While the SFRs for the [O iii] λ4363-detected and [O iii] λ4363-non-detected galaxies are modest (≈0.1–10 M_☉ yr⁻¹), their stellar masses are 1–2 dex lower than galaxies generally observed at $z\sim 1$. Therefore, we find that our emission-line galaxies are undergoing relatively strong star formation. The specific SFRs (SFR per unit stellar mass, SFR/M_⋆; hereafter sSFR) that we measure are between ${10}^{-10.8}$ and ${10}^{-6.1}$ yr⁻¹ with an average of 10${}^{-8.4}$ yr⁻¹ for the [O iii] λ4363-detected sample, and between ${10}^{-10.4}$ and ${10}^{-6.9}$ yr⁻¹ with an average of ${10}^{-8.8}$ yr⁻¹ for the [O iii] λ4363-non-detected sample. These averages are illustrated in Figure 2 by the dashed black line and dotted black line for the [O iii] λ4363-detected and [O iii] λ4363-non-detected samples, respectively. The gray shaded regions indicate the 1σ dispersion in sSFR for the samples.

Zoom In Zoom Out Reset image size

These sSFRs are enhanced by 0.25–4.0 dex above the M_⋆–SFR relation for $z\sim 0$ SDSS galaxies (Salim et al. 2007). Extrapolating the M_⋆–SFR relation of Whitaker et al. (2014a) and de los Reyes et al. (2015) toward lower stellar mass, we find that the sSFRs of our emission-line galaxies are ≈0.0–3.0 dex higher than "typical" galaxies at $z\sim 0.45$–0.85. While our sample is biased toward stronger star formation activity (see Section 4.1), 44% of [O iii] λ4363-detected and [O iii] λ4363-non-detected galaxies lie within ±0.3 dex (i.e., 1σ) of the $z\sim 0.8$ M_⋆–SFR relation of de los Reyes et al. (2015), and an additional 17% of our samples lie below the M_⋆–SFR relation by more than 0.3 dex (see Figure 2). For comparison, our previous [O iii] λ4363-detected study (Ly et al. 2014), which had shallower spectroscopy by a factor of ∼2, yielded a significant sSFR offset of $\approx 1$ dex on the M_⋆–SFR relation from typical star-forming galaxies at $z\leqslant 1$. The deeper observations of ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ result in a lower sSFR by $\approx 0.5$ dex. We also illustrate the Ly et al. (2015) [O iii] λ4363-selected metal-poor sample from DEEP2 as blue squares and triangles in Figure 2. This DEEP2 sample consists of galaxies with higher SFR activity than our [O iii] λ4363-detected and [O iii] λ4363-non-detected samples, which is in part due to the shorter integration time of DEEP2 (1 hr) than ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ (2 hr). It can also be seen that the ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ sample extends to lower stellar mass by $\approx 1$ dex at $z\sim 1$ than Ly et al. (2015).

We illustrate in Figure 3 the R₂₃ and O₃₂ strong-line ratios (Pagel et al. 1979):

$$\begin{eqnarray}&&{R}_{23}\equiv \displaystyle \frac{[{\rm{O}}\,{\rm{II}}]\,\lambda \lambda 3726,3729+[{\rm{O}}\,{\rm{III}}]\,\lambda \lambda 4959,5007}{{\rm{H}}\beta },\ \text{and}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{O}_{32}\equiv \displaystyle \frac{[{\rm{O}}\,{\rm{III}}]\,\lambda \lambda 4959,5007}{[{\rm{O}}\,{\rm{II}}]\,\lambda \lambda 3726,3729}.\end{eqnarray} \tag{ 5 }$$

We compare our [O iii] λ4363-detected and [O iii] λ4363-non-detected samples to typical $z\sim 2$ star-forming galaxies identified by the "MOSDEF" survey (black diamonds in this figure; Kriek et al. 2015; Shapley et al. 2015). We find that our metal-poor galaxies have similar interstellar properties (low metallicity, high ionization parameter) to the higher redshift galaxy population, suggesting that we have identified low-z analogs to $z\gtrsim 2$ galaxies. Specifically, the MOSDEF survey detects [O iii]λ5007 at S/N = 100 of ≈3 × 10⁻¹⁶ erg s⁻¹ cm⁻² or a line luminosity of 4 × 10⁴² erg s⁻¹ at z = 1.5, $1.3\times {10}^{43}$ erg s⁻¹ at z = 2.35, and 3 × 10⁴³ erg s⁻¹ at z = 3.35 (Kriek et al. 2015). As illustrated in Figure 24 of Paper I, the average [O iii] luminosity of the [O iii] λ4363-detected sample from ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ is 1.3–2.1 dex lower than the sensitivity of MOSDEF. Because the MOSDEF survey integrated for ∼1–2 hr, the [O iii] λ4363 emission for galaxies at $z\gtrsim 1.3$ would require at least ∼100 hr of Keck/MOSFIRE observations for individual S/N = 3 detections.

Zoom In Zoom Out Reset image size

We illustrate in Figure 4 the dependence of oxygen abundance on stellar mass in three redshift bins, $z\leqslant 0.3$, z = 0.3–0.5, and z = 0.5–1. In Figure 5, we compare the [O iii] λ4363-detected and [O iii] λ4363-non-detected samples from this paper and the DEEP2 [O iii] λ4363-detected and [O iii] λ4363-non-detected samples from Ly et al. (2015) against the Andrews & Martini (2013; hereafter AM13) M_⋆–Z relation of the form:

$$\begin{eqnarray}&&12+\mathrm{log}({\rm{O/H}})=12+\mathrm{log}{({\rm{O/H}})}_{{\rm{asm}}}-\mathrm{log}\left[1+{\left(\displaystyle \frac{{M}_{{\rm{TO}}}}{{M}_{\star }}\right)}^{\gamma }\right],\\ \end{eqnarray} \tag{ 6 }$$

where 12 + $\mathrm{log}({\rm{O/H}})$_asm is the asymptotic metallicity at the high mass end, ${M}_{{\rm{TO}}}$ is the turnover mass or "knee" in the M_⋆–Z relation, and γ is the slope of the low-mass end. This formalism is consistent with Moustakas et al. (2011) in describing the M_⋆–Z relation, and provides an intuitive understanding for the shape of the M_⋆–Z relation. For $z\sim 0.1$, AM13 find a best fit of 12 + $\mathrm{log}({\rm{O/H}})$ ${}_{{\rm{asm}}}=8.798$, $\mathrm{log}({M}_{{\rm{TO}}}/{M}_{\odot })=8.901$, and $\gamma =0.640$. At a given stellar mass, these emission-line selected samples are (on average) offset in 12 + $\mathrm{log}({\rm{O/H}})$ by ${0.13}_{-0.07}^{+0.06}$ dex at $z\leqslant 0.3$, −0.17${}_{-0.03}^{+0.07}$ dex at z = 0.3–0.5, and −0.24 ± 0.03 dex at z = 0.5–1. This demonstrates a moderate evolution in the M_⋆–Z relation of:

$$\begin{eqnarray}&&12\,+\,\mathrm{log}{({\rm{O/H}})-Z({M}_{\star })}_{\mathrm{AM}13}=A+B\mathrm{log}(1+z),\end{eqnarray} \tag{ 7 }$$

where $A={0.29}_{-0.13}^{+0.04}$ and $B=-{2.32}_{-0.26}^{+0.52}$. To better understand this evolution, we compute the average and median in each stellar mass bin, provided in Table 1, and shown as brown squares (average) and circles (median) in Figure 4. We then fit the averages with Equation (6) using MPFIT (Markwardt 2009). The fitting is repeated 10,000 times with each fit using the bootstrap approach to compute the average in each stellar mass bin. The best-fitting results are provided in Table 2 and the confidence contours are illustrated in Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Binned M_⋆–Z Relations

$\mathrm{log}({M}_{\star }/{M}_{\odot })$	N	$\langle Z\rangle$	Median Z	${\sigma }_{{\rm{obs}}}$	${\sigma }_{{\rm{int}}}$
(dex)		(dex)	(dex)	(dex)	(dex)
(1)	(2)	(3)	(4)	(5)	(6)
$z\leqslant 0.3$ (${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ Only)
7.25 ± 0.25	8	${7.85}_{-0.12}^{+0.12}$	${7.83}_{-0.09}^{+0.34}$	0.35	0.30
7.75 ± 0.25	6	${8.16}_{-0.09}^{+0.07}$	${8.12}_{-0.06}^{+0.23}$	0.21	0.19
8.25 ± 0.25	3	${8.53}_{-0.11}^{+0.10}$	${8.46}_{-0.05}^{+0.00}$	0.17	0.17
8.75 ± 0.25	7	${8.26}_{-0.14}^{+0.12}$	${8.22}_{-0.05}^{+0.18}$	0.37	0.36
$0.3\lt z\leqslant 0.5$ (${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ Only)
7.25 ± 0.25	5	${7.90}_{-0.10}^{+0.13}$	${7.95}_{-0.11}^{+0.06}$	0.27	0.25
7.75 ± 0.25	6	${7.90}_{-0.15}^{+0.12}$	${8.02}_{-0.28}^{+0.22}$	0.36	0.30
8.25 ± 0.25	12	${8.14}_{-0.08}^{+0.09}$	${8.24}_{-0.14}^{+0.00}$	0.36	0.34
8.75 ± 0.25	19	${8.27}_{-0.07}^{+0.06}$	${8.30}_{-0.07}^{+0.05}$	0.30	0.29
9.25 ± 0.25	10	${8.15}_{-0.12}^{+0.13}$	${8.19}_{-0.19}^{+0.16}$	0.41	0.39
$0.5\lt z\leqslant 1.0$ (${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ + Ly et al. 2015)
7.50 ± 0.25	5	${7.63}_{-0.14}^{+0.11}$	${7.58}_{-0.20}^{+0.00}$	0.26	0.22
8.00 ± 0.25	19	${7.92}_{-0.07}^{+0.07}$	${8.04}_{-0.07}^{+0.03}$	0.32	0.29
8.50 ± 0.25	32	${8.10}_{-0.05}^{+0.04}$	${8.10}_{-0.06}^{+0.09}$	0.25	0.23
9.00 ± 0.25	37	${8.25}_{-0.04}^{+0.04}$	${8.26}_{-0.02}^{+0.04}$	0.25	0.22
9.50 ± 0.25	11	${8.41}_{-0.08}^{+0.09}$	${8.38}_{-0.03}^{+0.05}$	0.28	0.26
10.00 ± 0.25	3	${8.37}_{-0.16}^{+0.23}$	${8.49}_{-0.00}^{+0.16}$	0.38	0.38

Notes. (1): Stellar mass bin. (2): Number of galaxies in each stellar mass bin, N. (3): Average 12 + $\mathrm{log}({\rm{O/H}})$. (4): Median 12 + $\mathrm{log}({\rm{O/H}})$. (5): Observed dispersion in 12 + $\mathrm{log}({\rm{O/H}})$ in each stellar mass bin. (6): Intrinsic dispersion in 12 + $\mathrm{log}({\rm{O/H}})$ after accounting for the average 12 + $\mathrm{log}({\rm{O/H}})$ measurement uncertainty, ${\sigma }_{{\rm{int}}}=\sqrt{{\sigma }_{{\rm{obs}}}^{2}-{\langle {\rm{\Delta }}{\rm{O/H}}\rangle }^{2}}$. Uncertainties for averages and medians are reported at the 16th and 84th percentile. These uncertainties are determined by statistical bootstrapping: random sampling with replacement, repeated 10,000 times.

Download table as:  ASCIITypeset image

Table 2.  Best Fit to Binned M_⋆–Z Relations

Redshift	12 + $\mathrm{log}({\rm{O/H}})$asm	$\mathrm{log}({M}_{{\rm{TO}}}/{M}_{\odot })$	γ
	(dex)	(dex)
(1)	(2)	(3)	(4)
$z\sim 0.1$ a	8.798	8.901	0.64
$z\leqslant 0.3$	${8.78}_{-0.10}^{+0.11}$	8.901b	${0.47}_{-0.11}^{+0.11}$
$0.3\lt z\leqslant 0.5$	${8.49}_{-0.07}^{+0.07}$	8.901b	${0.29}_{-0.09}^{+0.10}$
$0.5\lt z\leqslant 1.0$	${8.53}_{-0.05}^{+0.06}$	8.901b	${0.57}_{-0.08}^{+0.10}$
$0.5\lt z\leqslant 1.0$	${8.46}_{-0.05}^{+0.34}$	${8.61}_{-0.60}^{+0.57}$	${0.67}_{-0.09}^{+0.30}$

Notes. (1): Redshift. (2): Asymptotic metallicity at the high stellar mass end of the M_⋆–Z relation. (3): Turnover mass in the M_⋆–Z relation. (4): Slope of the low-mass end of the M_⋆–Z relation. See Equation (6) and Section 3.4 for further information. Uncertainties are reported at the 16th and 84th percentile, and are determined from the probability functions marginalized over the other two fitting parameters. Figure 6 illustrates the confidence contours for all three fitting parameters.

^aFrom AM13. ^bThe turnover mass was fixed to the value from AM13.

Download table as:  ASCIITypeset image

With only four or five stellar mass bins below M_⋆ ∼ 10⁹ M_☉ for our two lowest redshift bins ($z\lt 0.5$), fitting results are poorly constrained with all three parameters free. Specifically, the turnover mass (${M}_{{\rm{TO}}}$) and the asymptotic metallicity (12 +$\mathrm{log}({\rm{O/H}})$_asm) require measurements at higher stellar masses. In addition, with smaller sample sizes for these redshifts the best fits can easily be affected by a small number of outliers (e.g., the highest mass bin for $0.3\lt z\lt 0.5$ has a large number of metal-poor galaxies). For these reasons, we fixed ${M}_{{\rm{TO}}}$ to the local value obtained by AM13: ${10}^{8.901}$ M_☉. We find that the best-fit result to the SDF M_⋆–Z relation at $z\lt 0.3$ is consistent with AM13 within measurement uncertainties. At $0.3\lt z\lt 0.5$, the best fit yields a lower 12 + $\mathrm{log}({\rm{O/H}})$_asm by ≈0.3 dex. At z = 0.5–1, our observations extend to M_⋆ ∼ 10¹⁰ M_☉ and there are significantly more galaxies to better constrain the shape. Thus, we allow ${M}_{{\rm{TO}}}$ to be a free parameter in the fitting, in addition to adopting the local value. Our best fits to the M_⋆–Z relation at z = 0.5–1 indicate that the shape of the M_⋆–Z relation remains unchanged at $z\sim 1$, but with a lower 12 + $\mathrm{log}({\rm{O/H}})$_asm by $\approx 0.30$ dex (i.e., a lower metallicity at all stellar masses).

Several observational and theoretical investigations have proposed that the M_⋆–Z relation has a secondary dependence on the SFR (see e.g., Ellison et al. 2008; Lara-López et al. 2010; Mannucci et al. 2010; Davé et al. 2011; Lilly et al. 2013; Salim et al. 2014, and references therein). Specifically, the lower abundances at higher redshift may be explained by higher sSFR, such that there is a non-evolving (i.e., fundamental) relation (Lara-López et al. 2010; Mannucci et al. 2010). To test this relation, we adopt a non-parametric method of projecting the M_⋆–Z–SFR relation in various two-dimensional spaces.

First, we illustrate in Figure 7 the location on the M_⋆–SFR plane (Noeske et al. 2007; Salim et al. 2007) for galaxies in five different metallicity bins. We then compute the average and median sSFR for each bin. These are shown as brown and black solid lines, respectively, in each panel, and are summarized in the lower right panel. The hypothesis we are testing is whether, for a given stellar mass, galaxies shift toward higher SFRs as metallicity decreases. Our results show that indeed the sSFR is lower for higher values of $\mathrm{log}({\rm{O/H}})$, except at the lowest abundance bin. For our lowest abundance bin, Figure 7 shows that the distribution in sSFR is skewed (as evident by a ∼0.2 dex difference between the median and average values) by a small number of higher mass galaxies with low sSFR. The $\mathrm{log}({\rm{sSFR}})$–$\mathrm{log}({\rm{O/H}})$ slope that we measure is shallower than AM13, –0.30 versus –1.80. Specifically, the greatest difference in sSFR of ∼0.8 dex is at high metallicities. This difference is likely caused by a bias in our survey toward higher sSFR because metal-rich galaxies with low SFRs will not have [O iii] λ4363 detections and will fall below our flux limit cuts adopted for the [O iii] λ4363-non-detected sample. We defer a discussion on selection bias to Section 4.1.

Zoom In Zoom Out Reset image size

Next, we consider a projection first adopted by Salim et al. (2014): O/H as a function of the vertical offset on the M_⋆–SFR relation. The offset, defined as ${\rm{\Delta }}{({\rm{sSFR}})}_{{\rm{MS}}}$, measures the excess of star formation relative to "normal" galaxies of the same stellar mass and redshift. To facilitate comparisons with the local results of AM13, we use the Salim et al. (2007) $z\sim 0.1$ M_⋆–SFR relation as our reference relation:

$$\begin{eqnarray}&&\mathrm{log}({\rm{SFR}}/{M}_{\odot }\,{{\rm{yr}}}^{-1})=0.65\mathrm{log}({M}_{\star }/{M}_{\odot })-6.33.\end{eqnarray} \tag{ 8 }$$

This metallicity–${\rm{\Delta }}{({\rm{sSFR}})}_{{\rm{MS}}}$ comparison is performed in different stellar mass bins, and is illustrated in Figure 8. Here, we compare our sample to AM13, which is indicated by filled gray squares. This figure illustrates that our emission-line galaxy samples are qualitatively consistent with AM13; however, the sSFR dependence is weak. Specifically, there is a shallow inverse dependence at intermediate stellar masses ($8.1\leqslant \mathrm{log}({M}_{\star }/{M}_{\odot })\lt 8.6$), but no significant dependence in the remaining stellar mass bins. For these other stellar mass bins, there may be evidence for a positive metallicity–${\rm{\Delta }}{({\rm{sSFR}})}_{{\rm{MS}}}$ dependence; however, this is weak, with a significant dispersion of $\approx 0.3$ dex that is larger than measurement uncertainties.

Zoom In Zoom Out Reset image size

The last projection that we consider is how the M_⋆–Z relation depends on sSFR. This is illustrated in Figure 9 in five different sSFR bins from $\mathrm{log}({\rm{sSFR}}/{{\rm{yr}}}^{-1})=-9.8$ to –6.4. Similar to Figures 7 and 8, we overlay the AM13 sample as filled gray squares. The lower right panel of Figure 9 illustrates the median (black points) and average (red points) metallicities relative to the AM13 M_⋆–Z relation (see Equation (6)). While we find good agreement with AM13 at $-9.00\lt {\rm{sSFR}}\lt -8.25$, our results are broadly inconsistent with theirs. Specifically, we find that the relative offset on the M_⋆–Z relation increases with increasing sSFR. However, as discussed earlier, our selection function misses metal-rich galaxies. This effect has the largest impact for the lowest sSFR bin; the upper left panel of Figure 9 shows that metal-rich galaxies at M_⋆ ≳ 10⁹ M_☉ are not included in our sample. The effect of this selection would shift the average and median metallicities lower, possibly producing a false positive dependence. Because of the inability to measure metallicity in metal-rich galaxies with low star formation activity, the use of spectral stacking (such as AM13) is necessary to obtain average T_e-based abundances.

Zoom In Zoom Out Reset image size

One concern with our spectroscopic survey is the selection bias of requiring the detection of [O iii] λ4363. Specifically, detection of this line primarily depends on the electron temperature (or gas metallicity), which corresponds to

$$\begin{eqnarray}&&{T}_{e}([{\rm{O}}\,{\rm{III}}])=a{(-\mathrm{log}({ \mathcal R })-b)}^{-c},\text{where}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{ \mathcal R }\equiv \displaystyle \frac{F([{\rm{O}}\,{\rm{III}}]\,\lambda 4363)}{F([{\rm{O}}\,{\rm{III}}]\,\lambda \lambda 4959,5007)},\end{eqnarray} \tag{ 10 }$$

a = 13205, b = 0.92506, and c = 0.98062 (Nicholls et al. 2014), and the dust-corrected SFR and redshift, which determine the emission-line fluxes. At high SFRs, the probability of detecting [O iii] λ4363 is greater for a wide range of metallicity. This range in metallicity reduces such that only metal-poor galaxies with low SFRs can be detected in an emission-line flux limited survey.

To assess the selection function of our study, we examine the detectability of [O iii] λ4363 with MMT and Keck as a function of redshift, metallicity, Hβ luminosity (i.e., SFR), and dust attenuation. To determine the ${ \mathcal R }$ line ratio, we adopt a relation between T_e and 12 + $\mathrm{log}({\rm{O/H}})$ that is empirically based on our sample of [O iii] λ4363 detections:

$$\begin{eqnarray}&&{t}_{3}=28.767-5.865x+0.306{x}^{2},\end{eqnarray} \tag{ 11 }$$

where ${t}_{3}\equiv {T}_{e}$([O iii])/10${}^{4}\,{\rm{K}}$ and $x\equiv 12\,+\,\mathrm{log}({\rm{O/H}})$. This O/H–T_e relation is similar to that of Nicholls et al. (2014), which is based on several local samples (see their Figure 2).⁶ Then we determine the [O iii] λ5007/Hβ flux ratio as a function of metallicity by adopting $\mathrm{log}({{\rm{O}}}^{+}/{{\rm{O}}}^{++})=-0.114$, and the following equation from Izotov et al. (2006):

$$\begin{eqnarray}&&12+\mathrm{log}\left(\displaystyle \frac{{{\rm{O}}}^{++}}{{{\rm{H}}}^{+}}\right)\\ &&=\mathrm{log}\left[\displaystyle \frac{F([{\rm{O}}\,{\rm{III}}]\,\lambda \lambda 4959,5007)}{F({\rm{H}}\beta )}\right]\\ &&\,+6.200+\displaystyle \frac{1.251}{{t}_{3}}-0.55\mathrm{log}{t}_{3}-0.014{t}_{3}.\end{eqnarray} \tag{ 12 }$$

This value of $\mathrm{log}({{\rm{O}}}^{+}/{{\rm{O}}}^{++})$ is the average of our [O iii] λ4363-detected sample, which does not appear to be dependent on T_e across 10⁴–$2.5\times {10}^{4}$ K. We also examine local galaxies from Berg et al. (2012) and find no evidence for a T_e dependence across 10⁴–2 × 10⁴ K with an average $\mathrm{log}({{\rm{O}}}^{+}/{{\rm{O}}}^{++})=-0.166$ that is similar to our measured average. The combination of ${ \mathcal R }$, [O iii] λ5007/Hβ flux ratio, Hβ luminosity, and redshift determines the [O iii] λ4363 line flux:

$$\begin{eqnarray}&&F([{\rm{O}}\,{\rm{III}}]\,\lambda 4363)={ \mathcal R }\displaystyle \frac{1.33F([{\rm{O}}\,{\rm{III}}]\,\lambda 5007)}{F({\rm{H}}\beta )}\displaystyle \frac{L({\rm{H}}\beta )}{4\pi {d}_{L}^{2}},\end{eqnarray} \tag{ 13 }$$

where d_L is the luminosity distance and the [O iii] λ5007/λ4959 flux ratio is $\approx 3$ (Storey & Zeippen 2000). We illustrate in Figure 10 the average 3σ [O iii] λ4363 line sensitivity for the MMT and Keck spectra. Here, the sensitivity is computed by measuring the rms in the continuum of the spectra. We illustrate the effects of dust attenuation on the expected [O iii] λ4363 line flux by considering three different E(B − V) values: 0.13 (the average in our [O iii] λ4363 sample), 0.0, and 0.26 (the range encompasses ±1σ). The curves of [O iii] λ4363 line sensitivity are computed from MMT (Keck) spectra for four (five) average redshifts and overlaid in this figure with different colors. Because the on-source exposure time varies by a factor of a few to several, we normalized all estimates to two hours of integration (t₀). The observed points in Figure 10 account for the individual integration times with an offset to the Hβ luminosity of $0.5\mathrm{log}({t}_{{\rm{int}}}/{t}_{0})$. Typically, it can be seen that our [O iii] λ4363-detected galaxies lie to the right of our line sensitivities, while the [O iii] λ4363-non-detected galaxies lie to the left of the line sensitivity.

Zoom In Zoom Out Reset image size

Figure 10 demonstrates that the sensitivity to detect [O iii] λ4363 at solar (half-solar) metallicity, assuming the same SFR or Hα or Hβ luminosity, is on average 2.9 (1.7) times lower than at 12 + $\mathrm{log}({\rm{O/H}})$ $\leqslant$ 8.0. This suggests that the high stellar mass end (above 10⁹ M_☉) of our M_⋆–Z relations is likely biased toward lower metallicity, and that the M_⋆–Z relations could be steeper than reported in Section 3.4. We note, however, that this bias is relatively modest (less than a factor of 3 in sensitivity); thus, stacking at least ∼25 MMT and Keck spectra of metal-rich galaxies will yield average detections of [O iii] λ4363 that are significant above S/N = 5 at $z\sim 1$. A forthcoming paper of ${ \mathcal M }{ \mathcal A }{ \mathcal C }{ \mathcal T }$ will explore measurements from stacking, and further examine our selection bias.

Narrowband-selected sample. One of the first studies to use the narrowband imaging technique to select high-EW emission-line galaxies to obtain T_e-based metallicity was Kakazu et al. (2007). They targeted narrowband-excess emitters in the GOODS fields with Keck/DEIMOS and obtained 23 galaxies with $\geqslant 3\sigma$ detection of [O iii] λ4363 (Hu et al. 2009). While our [O iii] λ4363-detected sample is similar to theirs, Hu et al. (2009) mostly measured T_e below 12 + $\mathrm{log}({\rm{O/H}})$ ∼ 8.0, whereas our sample spans a wider range in metallicity at a given stellar mass or M_B.

Magnitude-limited sample. Our SDF [O iii] λ4363 sample at z = 0.2–1 overlaps closely in the M_⋆–Z plane with those measured by DEEP2 at $z\sim 0.8$ (Ly et al. 2015, see Figure 4). The main apparent difference is that the DEEP2 galaxies are from a magnitude-limited (i.e., M_⋆-limited) sample, which selects galaxies above M_⋆ $\sim \,{10}^{8.5}$ M_☉, and therefore higher metallicity.

In contrast, the Amorín et al. (2014, 2015) [O iii] λ4363-detected samples from VUDS and zCOSMOS, respectively, are strongly biased to only metal-poor galaxies at all stellar masses. This is well-illustrated in Figure 4, where nearly all of their galaxies are below the median of our at galaxies z = 0.5–1 (solid brown line in the lower left panel). These surveys obtain spectra at a lower resolution ($R\sim 200$), which limits the sensitivity to the detection of weak emission lines, particularly detecting [O iii] λ4363 in more metal-rich galaxies. This strong selection explains why galaxies from the Amorín et al. (2014, 2015) samples have systematically lower metallicities than galaxies in other samples.

The Izotov et al. (2014) and Berg et al. (2012) samples at low redshift also show systematically lower O/H than our SDF galaxies at $z\lesssim 0.3$. Izotov et al. (2014) selected galaxies from SDSS with high Hβ EWs and at $z\sim 0.2$. Because SDSS is a shallow magnitude-limited survey, their sample selection biased them toward more massive galaxies (M_⋆ ≳ 10⁹ M_☉) with lower metallicities. Berg et al. (2012) also reported that their sample has lower abundances at a given stellar mass when compared with other local M_⋆–Z studies (Lee et al. 2006). It appears that our $z\lesssim 0.3$ M_⋆–Z relation is consistent with the M_⋆–Z relation of Lee et al. (2006); however, we find a steeper slope.

As discussed earlier, the shape and evolution of the M_⋆–Z relation are important constraints for galaxy formation models, because the heavy-element abundances are set by enrichment from star formation, with dilution and loss from gas inflows and outflows, respectively (see Somerville & Davé 2015, and references therein). Efforts have been made to predict the M_⋆–Z relation from large cosmological simulations that either hydrodynamically model the baryons in galaxies (e.g., Davé et al. 2011; Obreja et al. 2014; Vogelsberger et al. 2014; Ma et al. 2016; Schaye et al. 2015) or adopt semi-analytical models with simple prescriptions for the baryonic physics (Gonzalez-Perez et al. 2014; Lu et al. 2014; Porter et al. 2014; Henriques et al. 2015; Croton et al. 2016).

We examine how well these numerical models and simulations predict the M_⋆–Z relation in Figure 11. Here we compare $z\sim 0$ predictions against AM13 in the left panel, and compare $z\sim 1$ predictions against our sample of z = 0.5–1 galaxies in the right panel. We note that the normalization of these predictions for the M_⋆–Z relation is dependent on the nucleosythesis yield, which is only accurately measured to ∼50% (R. Davé & K. Finlator 2015, private communication). Thus, the normalization cannot be used to compare with observations. For this reason, we normalize all M_⋆–Z relation predictions to 12 + $\mathrm{log}({\rm{O/H}})$ = 8.5 at ${M}_{\star }={10}^{9}$ M_☉ (at $z\sim 0;$ consistent with AM13), and examine relative evolution. For simplicity and consistency, predictions from hydrodynamic models are indicated by the dashed lines while semi-analytical model predictions are denoted by the dotted–dashed lines.

Zoom In Zoom Out Reset image size

First, we consider the predictions from the vzw simulation by Davé et al. (2011), which adopts "momentum-conserving" stellar winds. Their result is illustrated in Figure 11 by the gray dashed lines with the gray shaded regions encompassing the 16th and 84th percentile. At M_⋆ $\sim {10}^{9}$ M_☉, the slope in their M_⋆–Z relation is consistent with what we and AM13 measure. However, there are two issues with their predictions: (1) the decline in abundances with redshift (from z = 0 to $z\sim 1$) that they measure (∼0.1 dex; see Figure 2 in Davé et al. 2011) is much lower than what we observe ($\approx 0.25$ dex). (2) They predict a steep M_⋆–Z relation at higher stellar masses at all redshifts. This was not seen in our observations or by AM13. Unfortunately, the models from Davé et al. (2011) are unable to probe galaxies below M_⋆ $\approx \,{10}^{8.4}$ M_☉, where a steepening of the M_⋆–Z slope is seen at $z\sim 0$ and $z\sim 1$.

Next, we compare our results against "zoom-in" hydrodynamical simulations from the MaGICC (Making Galaxies in a Cosmological Context; Brook et al. 2012) and FIRE (Feedback in Realistic Environments; Hopkins et al. 2014) projects. While these simulations consist of much fewer galaxies than Davé et al. (2011), they provide higher spatial resolution on individual galaxies to resolve the structure of the ISM, star formation, and feedback, and span a wider range in galaxy stellar masses. For MaGICC, the results from Obreja et al. (2014) are illustrated by the red–orange squares with the best linear fit shown by the red–orange dashed lines. For FIRE, the redshift-dependent linear function described in Ma et al. (2016) is illustrated by the red dashed lines in Figure 11 with red stars for individual galaxies.⁷ Relative to AM13, Obreja et al. (2014) measure a slightly steeper slope at M_⋆ ∼ 10⁹ M_☉ than that found by Davé et al. (2011), although this is within the uncertainties. At $z\sim 0.7$, Obreja et al. (2014) find that abundances are lower by $\approx 0.2$ dex than at $z\sim 0$, which is roughly consistent with our observed M_⋆–Z relation evolution. Similar to Davé et al. (2011), Ma et al. (2016) measures a slope that is consistent with our sample and the results from AM13 at M_⋆ ∼ 10⁹ M_☉. However, the FIRE simulations find that abundances are lower by $\approx 0.25$ dex at $z\sim 0.8$ than at $z\sim 0$, which is consistent with our results. Because of the limited sample size of the MaGICC and FIRE simulations (only 7 and 24 galaxies at $z\sim 1$, respectively), constraining the shape of the M_⋆–Z relation is difficult. Furthermore, the MaGICC (FIRE) simulations only have two (three) $z\sim 1$ dwarf galaxies below M_⋆ ∼ 10⁸ M_☉. These galaxies can provide the strongest constraints on stellar winds from the M_⋆–Z relation.

In addition, we also consider the predictions of Vogelsberger et al. (2014) from the Illustris simulations,⁸ which are overlaid in Figure 11 as the cyan dashed lines. Similar to Davé et al. (2011), Obreja et al. (2014), and Ma et al. (2016), they predict an M_⋆–Z slope that is consistent with AM13 and our $z\sim 1$ sample; however, much like Davé et al. (2011), the amount of evolution predicted at M_⋆ ∼ 10⁹ M_☉ is inconsistent with the $\approx 0.25$ dex evolution that we observe. Also, the steep M_⋆–Z slope at the high mass end that Vogelsberger et al. (2014) predict is inconsistent with our observations and those of AM13.

Finally, we also overlay the predictions from the EAGLE hydrodynamical simulations (Crain et al. 2015; Schaye et al. 2015) in Figure 11 as green dashed lines and green shaded regions encompassing the 16th and 84th percentile. The predictions are obtained from the public catalog (McAlpine et al. 2016).⁹ The EAGLE simulation results agree with the observed shape of the M_⋆–Z relation at ${M}_{\star }\gtrsim 3\times {10}^{8}$ M_☉ for $z\sim 0$ (AM13) and at ${M}_{\star }\gtrsim 3\times {10}^{8}$ M_☉ for $z\sim 1$. At lower stellar mass it predicts a shallow M_⋆–Z slope, which is inconsistent with observational results at $z\sim 0$. However, this shallow slope is believed to be caused by poor resolution because the turnover occurs when the number of star particles falls below ∼ 10⁴ (Schaye et al. 2015). The EAGLE simulation does predict ≈0.2 dex evolution in the M_⋆–Z normalization, which is consistent with our observed evolution.

For completeness, we also overlay the predictions from several semi-analytical models. These models adopt different assumptions and prescriptions. The predictions¹⁰ for Croton et al. (2016) are shown by the dotted–dashed orange line and orange diamonds¹¹ for $z\sim 0$ and the dotted–dashed orange line and orange shaded region for $z\sim 1$. We also overlay Gonzalez-Perez et al. (2014) as a black dotted–dashed line, Henriques et al. (2015) as the olive dotted–dashed line with olive shaded region indicating the 16th and 84th percentile, Lu et al. (2014) as the purple dotted–dashed line with purple shaded region indicating 68% dispersion, and Porter et al. (2014) as the yellow dotted–dashed line with black outlines.

We note that none of these semi-analytical models predict a moderate (≈0.25 dex) evolution in the M_⋆–Z relation at $z\lesssim 1$. They either find no evolution or no more than 0.1 dex. Croton et al. (2016), Henriques et al. (2015), and Porter et al. (2014) do predict a M_⋆–Z slope at M_⋆ ∼ 10⁹ M_☉ that is consistent with observations at $z\sim 0$ and $z\sim 1$. The semi-analytical model that disagrees significantly from observations is Lu et al. (2014). They predict the steepest M_⋆–Z slope and higher metallicities at a given stellar mass for $z\sim 1$. Croton et al. (2016) is able to reproduce the shape of the M_⋆–Z relation at $z\sim 0;$ however this is not a surprise because the local M_⋆–Z relation (Tremonti et al. 2004) is used as a secondary constraint in their model.

Given these comparisons, we find that the only models that can reproduce both the observed evolution in the M_⋆–Z relation ($\approx 0.25$ dex) and the slope at M_⋆ $\sim \,{10}^{9}$ M_☉ are the FIRE and EAGLE simulations. As discussed above, the FIRE simulations provide the highest spatial resolution to resolve stellar feedback (Hopkins et al. 2014), and the EAGLE simulation with the best agreement with the local M_⋆–Z relation has the highest particle resolution. While these results suggest that resolving the physical processes in the ISM is critical for further understanding the chemical enrichment process in galaxies, there has been some success in yielding the same M_⋆–Z relation at multiple epochs with different resolutions (Torrey et al. 2014). This would suggest that what may be more important for lower resolution simulations is correctly handling the baryonic gas flows on subresolution scales.

One reason why previous numerical simulations have not used the M_⋆–Z relation as an observational constraint is the growing concern that "strong-line" metallicity diagnostics may not be valid for higher redshift galaxies (e.g., Steidel et al. 2014; Sanders et al. 2015; Cowie et al. 2016; Dopita et al. 2016). However, now that the evolution of the temperature-based M_⋆–Z relation has been measured, we encourage forthcoming models and simulations to utilize the evolution of the T_e-based M_⋆–Z relation as an important constraint for galaxy formation models. In addition, we encourage future work to probe galaxies below M_⋆$\,\sim \,{10}^{8}$ M_☉, because this remains an unexplored parameter space where observations find a steep M_⋆–Z dependence (see Figure 11). Finally, should computing infrastructures allow, improving the particle resolution of large-scale galaxy simulations and using the "zoom-in" technique for more detailed studies are additional improvements that may provide a better understanding of the gas flows in galaxies.