The MOSDEF Survey: A Stellar Mass–SFR–Metallicity Relation Exists at z ∼ 2.3^∗

The MOSDEF survey was a 4-year program in which we utilized the Multi-Object Spectrometer For Infra-Red Exploration (MOSFIRE; McLean et al. 2012) on the 10 m Keck I telescope to obtain near-infrared (rest-frame optical) spectra of galaxies at 1.4 ≤ z ≤ 3.8 (Kriek et al. 2015). Galaxies were targeted in three redshift ranges such that strong optical emission lines fall in near-infrared windows of atmospheric transmission: 1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.61, and 2.95 ≤z ≤ 3.80. In this study, we focus on the middle redshift range at z ∼ 2.3. Targets were selected from the photometric catalogs of the 3D-HST survey (Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016) based on their redshifts and observed H-band (rest-frame optical) magnitudes. Spectroscopic or HST grism redshifts were used when available; otherwise photometric redshifts were utilized for selection. For the z ∼ 2.3 redshift bin, galaxies were targeted down to a fixed H-band AB magnitude of H = 24.5, as measured from HST/WFC3 F160W imaging. This targeting scheme effectively selects galaxies down to a rough stellar mass limit of log(M_*/M_⊙) ∼ 9.0, although the sample is not complete below log(M_*/M_⊙) ∼ 9.5 (see Shivaei et al. 2015). At 2.09 ≤ z ≤ 2.61, [O ii] λλ3726,3729 falls in the J band; Hβ and [O iii] λλ4959,5007 fall in the H-band; and Hα, [N ii] λλ6548,6584, and [S ii] λλ6716,6731 fall in the K band. Observations were completed in 2016 May, with the full survey having targeted ∼1500 galaxies and measured robust redshifts for ∼1300 galaxies, with ∼700 at z ∼ 2.3 and ∼300 at z ∼ 1.5 and 3.4, respectively. Full technical details of the MOSDEF survey observing strategy and data reduction can be found in Kriek et al. (2015).

2.2.1. Stellar Mass

2.2.2. Emission-line Fluxes and Redshift

2.2.3. Reddening Correction and Star Formation Rate

2.3.1. MOSDEF z ∼ 2.3 Sample

We selected a sample of star-forming galaxies from the full MOSDEF data set, with which we can investigate the existence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3. We required a secure redshift determination at 2.0 ≤ z ≤ 2.7, log(M_*/M_⊙) ≥ 9.0, and a detection of both Hα and Hβ at signal-to-noise ratio S/N ≥ 3 so that robust estimates of the reddening-corrected SFR may be obtained. The stellar mass cut removed six objects with lower stellar masses than the MOSDEF targeting scheme was designed to select. Four of these objects were not main targets but were instead serendipitously detected. Both the number of objects targeted and spectroscopic success rate drops precipitously below log(M_*/M_⊙) = 9.0 (Kriek et al. 2015). We rejected objects identified as AGN by their X-ray or infrared properties (Coil et al. 2015; Azadi et al. 2017), and additionally removed objects with [N ii] λ6584/Hα > 0.5 when [N ii] λ6584 is detected at S/N ≥ 3, which are not likely to be dominated by star formation. We additionally rejected objects with significant sky-line contamination of Hα or Hβ such that the line flux is unreliable. The requirement that Hβ is cleanly detected at S/N ≥ 3 was the most restrictive of these cuts, removing 45% of 462 galaxies that satisfy all other criteria. We discuss potential biases introduced by the Hβ detection requirement below. This selection yielded a sample of 260 galaxies at z_med = 2.29 with stellar masses spanning the range log(M_*/M_⊙) = 9.0–11.4 and SFRs ranging from 1.4 to 260 M_⊙ yr⁻¹. The median stellar mass and SFR of the sample is log(M_*/M_⊙) = 9.92 and 22 M_⊙ yr⁻¹, respectively.

The redshift distribution, and the SFR and specific SFR (sSFR = SFR/M_*) as a function of stellar mass, of the z ∼ 2.3 sample are shown in Figure 1. Galaxies in this sample scatter about the mean SFR–${M}_{* }$ relation at z ∼ 2.0 − 2.5, with the mean properties falling on the relations seen in past studies (e.g., Whitaker et al. 2012; Speagle et al. 2014; Shivaei et al. 2015). At log(M_*/M_⊙) < 9.5, galaxies below the mean SFR–${M}_{* }$ relation begin to fall below the MOSDEF Hβ detection threshold at z ∼ 2.3 (Kriek et al. 2015). Shivaei et al. (2015) showed that excluding objects with Hβ non-detections did not significantly bias a smaller sample of MOSDEF star-forming galaxies at log(M_*/M_⊙) > 9.5. Galaxies for which Hα is detected but Hβ is not detected, shown as SFR and sSFR lower limits in Figure 1 due to a lack of constraints on reddening, are distributed relatively uniformly with ${M}_{* }$. Above log(M_*/M_⊙) = 9.5, most limits lie well above the Hβ detection threshold at z ∼ 2.3, suggesting that most of the Hβ non-detections are due to sky-line contamination. Since sky-line contamination is essentially a uniform redshift selection that does not correlate with other galaxy properties, excluding these Hβ non-detections does not bias our sample (see Section 3.2 for additional discussion of Hβ selection effects). Shivaei et al. (2015) also demonstrated that below log(M_*/M_⊙) = 9.5, MOSDEF samples may be incomplete due to a bias against young objects with small Balmer and 4000 Å breaks for which the photometric redshifts used in our target selection may be inaccurate. This incompleteness does not significantly affect our results because the majority of our sample lies above log(M_*/M_⊙) = 9.5. The results of our analysis do not significantly change if we restrict the sample to log(M_*/M_⊙) > 9.5.

Zoom In Zoom Out Reset image size

We also investigated whether spectroscopic incompleteness biases our sample. Of ∼700 galaxies targeted at $2.0\leqslant z\leqslant 2.7$ in the MOSDEF survey, 93 failed to produce spectroscopic redshifts. Half of these (45/93) are quiescent galaxies based on their rest-frame UVJ colors and thus do not affect our sample of star-forming galaxies. Red star-forming galaxies at z ∼ 2.3 failed to yield redshifts 15% of the time, while redshifts were not measured for only 6% of blue star-forming galaxies. These statistics agree with what Kriek et al. (2015) reported based on one-third of the full MOSDEF data set. The spectroscopic failures do not occupy any region of UVJ color space in which there are no objects with measured redshifts. Our sample is thus slightly biased against red star-forming galaxies. However, our sample is dominated by blue star-forming galaxies (>85%) such that the bias against red star-forming galaxies has minimal impact on our results.

2.3.2. SDSS z ∼ 0 Comparison Sample

The MOSDEF data set crucially provides coverage of several strong rest-optical emission features, allowing us to measure multiple metallicity indicators widely applied in the local universe. We utilize multiple emission-line ratios to estimate metallicity in this analysis, including N2, O3N2, N2O2, and O32. We do not utilize R23 to estimate metallicity for reasons described below, but retain the use of this metallicity-sensitive line ratio empirically, along with O3.

For N2 and O3N2, we use the empirically calibrated relations of Pettini and Pagel (2004) based on individual z ∼ 0 H ii regions:

$\begin{eqnarray}&&12+\mathrm{log}({\rm{O/H}})=8.90+0.57\times {\rm{log(N2)}},\end{eqnarray} \tag{ 7 }$

and

$\begin{eqnarray}&&12+\mathrm{log}({\rm{O/H}})=8.73-0.32\times {\rm{log(O3N2)}}.\end{eqnarray} \tag{ 8 }$

The intrinsic uncertainties associated with these calibrations are 0.14 and 0.18 dex for N2 and O3N2, respectively. The N2 calibration is valid over the range −2.5 < log(N2) < −0.3, while the O3N2 calibration is applicable for log(O3N2) ≲ 2.0. All of the z ∼ 0 sample and the entire high-redshift MOSDEF sample except for six galaxies fall within these bounds.

Brown et al. (2016) provided a calibration of galaxy metallicity as a function of N2O2 and offset from the mean z ∼ 0 sSFR–${M}_{* }$ sequence (Δlog(sSFR)):

$\begin{eqnarray}12+\mathrm{log}({\rm{O/H}}) & = & 9.20+0.54\times {\rm{log(N2O2)}}\\ & & -0.36\times {\rm{\Delta }}{\rm{log(sSFR)}}.\end{eqnarray} \tag{ 9 }$

The Δlog(sSFR) term was included because these authors found that N2O2 systematically increases at fixed direct-method metallicity as Δlog(sSFR) increases. Sanders et al. (2017) demonstrated that the separation of galaxies of different Δlog(sSFR) in the N2O2–metallicity plane is a result of contamination of the integrated emission-line fluxes from diffuse ionized gas (DIG), which more strongly affects galaxies with lower Δlog(sSFR). After correcting the Brown et al. (2016) strong-line ratios and direct-method metallicities for DIG contamination, the dependence on Δlog(sSFR) disappears and the galaxies are shifted onto the H ii region metallicity scale. To estimate metallicities from N2O2, we perform a linear fit to the corrected Brown et al. (2016) data points from Sanders et al. (2017), obtaining the following expression:

$\begin{eqnarray}&&12+\mathrm{log}({\rm{O/H}})=8.94+0.73\times {\rm{log(N2O2)}}.\end{eqnarray} \tag{ 10 }$

The intrinsic scatter around this relation is 0.2 dex, and the calibration is valid over the range −1.3 < log(N2O2) < 0.0. The entire z ∼ 0 comparison sample falls in this range. Only 10 galaxies in the z ∼ 2.3 sample have log(N2O2) < −1.3, while the remainder fall in the range over which the calibration is valid.

We utilize the O32 calibration of Jones et al. (2015):

$\begin{eqnarray}&&12+\mathrm{log}({\rm{O/H}})=8.3439-0.4640\times {\rm{log(O32)}}.\end{eqnarray} \tag{ 11 }$

The intrinsic calibration uncertainty is 0.11 dex. This relation was calibrated using a sample of 113 z ∼ 0 star-forming galaxies with metallicity estimates based on electron temperatures, and no evolution in this calibration is seen out to z ∼ 0.8 (Jones et al. 2015). The calibrating data set only covers the range 0.0 ≲ log(O32) ≲ 1.0, with the majority of the data set lying at 0.0 < log(O32) < 0.5. While 91 galaxies in the z ∼ 2.3 sample and all of the local comparison sample has log(O32) < 0.0, other empirical and theoretical calibrations do not show significant changes to the slope of O32 calibrations below log(O32) = 0.0 (Maiolino et al. 2008; Curti et al. 2017). We therefore apply the Jones et al. (2015) O32 calibration to galaxies with −1.0 < log(O32) < 1.0 while noting that a systematic bias in metallicity estimates from this calibration may be introduced at log(O32) < 0.0.

While R23 is commonly employed as a metallicity indicator in the local universe, we do not estimate metallicity from R23. R23 is double-valued as a function of metallicity, with the turnover occurring at 8.0 ≲ 12 + log(O/H) ≲ 8.4 and 0.7 ≲ log(R23) ≲ 1.0 (e.g., McGaugh 1991; Kewley & Dopita 2002; Pilyugin & Thuan 2005). The majority of our z ∼ 2.3 sample lies in this region of parameter space, where small observational uncertainties lead to large uncertainties in metallicity and assigning an object to the upper or lower R23 branch is non-trivial. Breaking the upper and lower branch degeneracy is usually achieved using N2 or N2O2 in combination with R23 (e.g., Kewley & Dopita 2002). Star-forming galaxies at z ∼ 2.3 display significant offsets from the z ∼ 0 population in strong-line ratio diagrams involving lines of nitrogen (Masters et al. 2014; Steidel et al. 2014; Shapley et al. 2015; Sanders et al. 2016b), however, suggesting that the criterion for upper or lower R23 branch assignment based on N2 or N2O2 evolves with redshift. Due to these uncertainties in R23-based metallicities for the high-redshift sample, we do not convert R23 values to metallicities but instead utilize R23 in an empirical sense only.

While all 260 individual galaxies in the z ∼ 2.3 sample have detections of Hα and Hβ, some do not have S/N ≥ 3 for [N ii] λ6584, [O ii] λλ3726,3729, and/or [O iii] λ5007. The sample contains 143 N2 detections, 126 O3N2 detections, 118 N2O2 detections, 169 O32 and R23 detections, and 223 O3 detections. Each of the detected line-ratio subsets has median stellar mass in the range 9.9 < log(M_*/M_⊙) < 10.1 and median SFR within the range 22 < SFR/M_⊙ yr⁻¹ < 32. We incorporate information from galaxies with metal line non-detections using a stacking analysis as described below.

As mentioned previously, DIG is a significant contaminant of global galaxy spectra at z ∼ 0 and introduces systematic biases in galaxy emission-line ratios and metallicity estimates (Sanders et al. 2017). In order to obtain accurate metallicities from strong-line calibrations, both the measured line ratios and the metallicity calibrations must be free of the effects of DIG contamination. Relations calibrated using individual H ii  regions such as the Pettini and Pagel (2004) N2 and O3N2 calibrations are free from this issue. We have explicitly corrected the z ∼ 0 sample of Brown et al. (2016) to the H ii  region metallicity scale. While the calibrating data set of the Jones et al. (2015) O32 relation has not been corrected, we do not expect a high level of DIG contamination in this case. The galaxies selected by Jones et al. (2015) are highly star-forming and have high-sSFR and large Hα surface brightness. Oey et al. (2007) found that the importance of DIG to total galaxy emission decreases with increasing Hα surface brightness, such that DIG emission is negligible in starburst galaxies. Using the calibration from Sanders et al. (2017) based on the Oey et al. (2007) data set, we estimate that the fraction of total Balmer emission originating from DIG (${f}_{\mathrm{DIG}}$) is ≲30% for the galaxies in the Jones et al. (2015) sample based on their Hα surface brightnesses. This result suggests that H ii region emission dominates the Jones et al. (2015) calibrating sample and the correction for DIG contamination would be minor compared to what is required for typical z ∼ 0 galaxies with ${f}_{\mathrm{DIG}}$ ≈ 55% (Sanders et al. 2017).

We corrected the z ∼ 0 comparison sample for DIG contamination using the results of Sanders et al. (2017), which are based upon empirical models treating galaxies as ensembles of H ii and DIG regions. These models also correct for flux-weighting biases present in integrated galaxy measurements because the line-emitting regions falling in the spectroscopic aperture have a range of metallicites. The AM13 ${M}_{* }$ stacks were corrected using the z ∼ 0 model with ${f}_{\mathrm{DIG}}$ = 0.55. For the AM13 SFR–${M}_{* }$ bins, Sanders et al. (2017) determined the median Hα surface brightness in each bin and constructed a model matched to the corresponding f_DIG, since the relative importance of DIG varies with SFR. These corrections allow us to obtain more robust metallicity estimates for the z ∼ 0 comparison sample by using the corrected line ratios with the metallicity calibrations described above.

There are currently no observational constraints on the importance of DIG at high redshifts. However, if high-redshift galaxies follow a similar relation between f_DIG and Hα surface brightness as local galaxies, the high SFRs and compact sizes of high-redshift galaxies suggest that DIG emission is negligible at z ∼ 2.3. We operate under this assumption, although high-resolution spatially resolved emission-line observations of high-redshift galaxies are ultimately required to reveal the role of DIG in the early universe. In addition to compact sizes and high SFRs, galaxies have much larger gas fractions at fixed ${M}_{* }$ at z ∼ 2 than at z ∼ 0 (Tacconi et al. 2013), and likely display differences in the distribution and importance of different ISM phases. In light of such differences, the assumption that the relation between f_DIG and Hα surface brightness does not change with redshift may need to be revisited.

We selected a subset of the sample of 260 z ∼ 2.3 star-forming galaxies with which we perform a stacking analysis. For this stacking subsample, we additionally required that the spectrum of each object had wavelength coverage of [O ii] λλ3726,3729, Hβ, [O iii] λ5007, Hα, and [N ii] λ6584, and that the stellar mass fell in the range 9.0 ≤ log(M_*/M_⊙) ≤ 10.5. The wavelength coverage criterion allows for the measurement of all emission-line ratios of interest in this study from the stacked spectra. The stellar mass criterion is designed to only stack over the mass range in which the MOSDEF z ∼ 2.3 star-forming sample is not obviously incomplete based on target selection and/or spectroscopic success rate. Kriek et al. (2015) showed that both the number of galaxies targeted and the spectroscopic success rate sharply decrease at log(M_*/M_⊙) < 9.0 using data from the first 2 years of the MOSDEF survey. Such behavior is expected because of the target selection down to fixed H-band magnitude and the inherent faintness of low-mass targets. A drop in the spectroscopic success rate at log(M_*/M_⊙) > 10.5 was also identified and found to be partially caused by a significantly lower success rate for red star-forming and quiescent galaxies compared to blue star-forming galaxies (Kriek et al. 2015). We have confirmed these trends in success rate with the full MOSDEF z ∼ 2.3 sample. Applying the above criteria yields a stacking subsample of 195 galaxies at z_med = 2.29 with median stellar mass log(M_*/M_⊙) = 9.89 and median SFR of 18 M_⊙ yr⁻¹. The stacking subsample has slightly lower median stellar mass and SFR than the larger z ∼ 2.3 sample it was selected from, primarily because of the upper stellar mass cut.

We used two different binning methods for this analysis. We binned in stellar mass in order to determine mean z ∼ 2.3 line-ratio and metallicity relations as a function of stellar mass. Galaxies were stacked in four stellar mass bins selected such that an equal number of galaxies fell in each bin (48–49 galaxies per bin). We refer to the stacks binned in stellar mass only as the "${M}_{* }$ stacks." We also binned in both stellar mass and offset from the mean z ∼ 2.3 sSFR–${M}_{* }$ relation (Δlog(sSFR)) in order to assess the presence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 and its evolution from z ∼ 0. We found that dividing the stacking sample at Δlog(sSFR/yr⁻¹) = ±0.2 splits the sample into three roughly equal parts, with 68 galaxies in the highly star-forming subsample with Δlog(sSFR) > 0.2, 67 galaxies falling within 0.2 dex of the mean relation, and 60 galaxies with Δlog(sSFR) < −0.2. Since these bins were selected with respect to the mean z ∼ 2.3 sSFR–${M}_{* }$ relation, the stellar mass distribution in each bin is similar except for the Δlog(sSFR) < −0.2 bin, in which objects with log(M_*/M_⊙) ≲ 9.5 fall below the MOSDEF Hβ detection limit at z ∼ 2.3. The Δlog(sSFR) < −0.2 bin is thus biased toward higher ${M}_{* }$. After dividing the sample in Δlog(sSFR), galaxies in each Δlog(sSFR) bin were stacked in two stellar mass bins divided at the median stellar mass, yielding a total of 6 ${M}_{* }$–Δlog(sSFR) bins. We refer to the stacks binned in both ${M}_{* }$ and Δlog(sSFR) as the "${M}_{* }$–ΔsSFR stacks."

We produced composite spectra by first shifting each spectrum into the rest frame and converting from flux density to luminosity density using the spectroscopic redshift. Since some of the line ratios of interest require reddening correction, we corrected each individual spectrum for reddening by applying a reddening correction to the luminosity density at each wavelength element based on the Balmer decrement (Hα/Hβ) of each galaxy, assuming the Cardelli et al. (1989) extinction curve. This process effectively only corrects the nebular lines, since the continuum is not significantly detected for individual galaxies in our stacking sample. Each reddening-corrected spectrum was then normalized by the dust-corrected Hα luminosity. Spectra were interpolated onto a wavelength grid with spacing equal to the wavelength sampling at the median redshift of the sample, equal to 0.40 Å in J band, 0.50 Å in H band, and 0.66 Å in K band. Normalized composite spectra were created by taking the median value of the normalized spectra at each wavelength increment, and multiplying by the median dust-corrected Hα luminosity to produce the final composite spectrum in units of luminosity density.

Emission-line luminosities were measured from stacked spectra by fitting a flat continuum to regions around emission features and Gaussian profiles to the emission lines following the fitting methodology for individual spectra. Uncertainties on emission-line fluxes were estimated using a Monte Carlo technique in which we perturbed the stellar masses according to the uncertainties and divided the perturbed sample into the same stellar mass bins as before, bootstrap resampled each bin to account for sample variance, perturbed the individual spectra according to the error spectra, perturbed the individual E(B–V)_gas values and applied them to correct the individual spectra for reddening, stacked according to the method described above, and remeasured the emission-line luminosities. This process was repeated 100 times, and the uncertainties on line luminosities were obtained from the 68th-percentile width of the distribution of remeasured line luminosities. The median Balmer absorption luminosity of the individual galaxies in each bin was utilized to apply Balmer absorption corrections to the stacked Hα and Hβ luminosities.

Emission-line ratios and uncertainties were calculated using the line luminosities and uncertainties measured from the stacked spectra. We found that the stacking method described above robustly reproduces line ratios characteristic of the individual galaxies in a stack by comparing the stacked line ratios to median line ratios of individual galaxies with significant detections of all of the emission lines of interest. Stellar mass and SFR were assigned to each stack using the median ${M}_{* }$ and median SFR of the individual galaxies in each bin. Uncertainties on the median ${M}_{* }$ and SFR were estimated using the distribution of median values produced by the Monte Carlo simulations. The median ${M}_{* }$ and SFR of galaxies in each bin and measured line ratios of the ${M}_{* }$ and ${M}_{* }$–ΔsSFR stacks are presented in Table 1.

Table 1.  Galaxy Properties and Emission-line Ratios from Stacks of z ∼ 2.3 Star-forming Galaxy Spectra

Ngala	$\mathrm{log}\left(\tfrac{{M}_{* }}{{M}_{\odot }}\right)$ b	SFRmedc	$\mathrm{log}\left(\tfrac{\mathrm{sSFR}}{{{yr}}^{-1}}\right)$ d	log(N2)	log(O3)	log(O3N2)	log(N2O2)	log(O32)	log(R23)
M* stacks
49	9.0−9.62; ${9.43}_{-0.03}^{+0.01}$	${12.1}_{-1.2}^{+3.4}$	$-{8.34}_{-0.05}^{+0.16}$	$-{1.22}_{-0.05}^{+0.10}$	${0.66}_{-0.02}^{+0.04}$	${1.89}_{-0.10}^{+0.06}$	$-{1.13}_{-0.08}^{+0.20}$	${0.28}_{-0.13}^{+0.10}$	${0.93}_{-0.03}^{+0.04}$
49	9.62−9.89; ${9.78}_{-0.01}^{+0.02}$	${10.6}_{-4.4}^{+0.8}$	$-{8.75}_{-0.15}^{+0.03}$	$-{1.08}_{-0.09}^{+0.06}$	${0.55}_{-0.02}^{+0.05}$	${1.63}_{-0.07}^{+0.11}$	$-{1.12}_{-0.18}^{+0.09}$	${0.07}_{-0.11}^{+0.16}$	${0.89}_{-0.05}^{+0.06}$
49	9.89−10.13; ${9.97}_{-0.04}^{+0.01}$	${16.6}_{-6.8}^{+2.8}$	$-{8.75}_{-0.13}^{+0.11}$	$-{0.90}_{-0.06}^{+0.05}$	${0.45}_{-0.04}^{+0.02}$	${1.36}_{-0.06}^{+0.06}$	$-{0.93}_{-0.14}^{+0.13}$	$-{0.05}_{-0.17}^{+0.07}$	${0.84}_{-0.05}^{+0.05}$
48	10.14−10.50; ${10.30}_{-0.02}^{+0.01}$	${45.1}_{-9.3}^{+14.9}$	$-{8.65}_{-0.07}^{+0.17}$	$-{0.71}_{-0.03}^{+0.03}$	${0.37}_{-0.05}^{+0.02}$	${1.09}_{-0.06}^{+0.03}$	$-{0.80}_{-0.15}^{+0.21}$	$-{0.20}_{-0.19}^{+0.23}$	${0.83}_{-0.16}^{+0.07}$
M*–ΔsSFR stacks: Δlog(sSFR/yr−1) < −0.2
30	9.18−9.90; ${9.72}_{-0.01}^{+0.04}$	${5.7}_{-1.5}^{+0.2}$	$-{8.96}_{-0.12}^{+0.01}$	$-{1.05}_{-0.11}^{+0.10}$	${0.43}_{-0.04}^{+0.04}$	${1.49}_{-0.14}^{+0.13}$	$-{1.04}_{-0.16}^{+0.13}$	${0.00}_{-0.11}^{+0.08}$	${0.80}_{-0.04}^{+0.05}$
30	9.90−10.49; ${10.04}_{-0.05}^{+0.02}$	${11.1}_{-0.6}^{+1.9}$	$-{8.99}_{-0.01}^{+0.10}$	$-{0.79}_{-0.06}^{+0.09}$	${0.34}_{-0.05}^{+0.04}$	${1.14}_{-0.10}^{+0.06}$	$-{0.75}_{-0.09}^{+0.13}$	$-{0.11}_{-0.12}^{+0.04}$	${0.76}_{-0.04}^{+0.05}$
M*–ΔsSFR stacks: $-0.2\leqslant {\rm{\Delta }}\mathrm{log}$(sSFR/yr ${}^{-1})\leqslant +0.2$
34	9.00−9.86;${9.50}_{-0.03}^{+0.01}$	${9.2}_{-2.1}^{+0.4}$	$-{8.53}_{-0.07}^{+0.03}$	$-{1.16}_{-0.17}^{+0.04}$	${0.59}_{-0.04}^{+0.03}$	${1.76}_{-0.07}^{+0.14}$	$-{1.09}_{-0.19}^{+0.10}$	${0.21}_{-0.13}^{+0.09}$	${0.88}_{-0.05}^{+0.03}$
33	9.86−10.47; ${10.1}_{-0.01}^{+0.05}$	${26.2}_{-0.1}^{+4.4}$	$-{8.68}_{-0.03}^{+0.07}$	$-{0.81}_{-0.04}^{+0.07}$	${0.41}_{-0.04}^{+0.04}$	${1.23}_{-0.08}^{+0.05}$	$-{0.83}_{-0.08}^{+0.11}$	$-{0.11}_{-0.09}^{+0.04}$	${0.84}_{-0.02}^{+0.05}$
M*–ΔsSFR stacks: Δlog(sSFR/yr−1) >+0.2
34	9.20−9.92; ${9.55}_{-0.05}^{+0.01}$	${32.1}_{-4.1}^{+3.0}$	$-{8.04}_{-0.04}^{+0.07}$	$-{1.23}_{-0.10}^{+0.07}$	${0.69}_{-0.01}^{+0.03}$	${1.93}_{-0.08}^{+0.11}$	$-{1.22}_{-0.11}^{+0.11}$	${0.25}_{-0.06}^{+0.15}$	${0.97}_{-0.04}^{+0.03}$
34	9.92−10.5; ${10.25}_{-0.03}^{+0.01}$	${88.2}_{-1.0}^{+17.3}$	$-{8.30}_{-0.02}^{+0.12}$	$-{0.78}_{-0.05}^{+0.05}$	${0.46}_{-0.03}^{+0.03}$	${1.24}_{-0.05}^{+0.08}$	$-{0.97}_{-0.07}^{+0.07}$	$-{0.18}_{-0.06}^{+0.09}$	${0.92}_{-0.04}^{+0.03}$

Notes.

^aNumber of galaxies in a bin. ^bRange and median $\mathrm{log}({M}_{* }/{M}_{\odot })$ of galaxies in a bin. ^cMedian dust-corrected Hα SFR of galaxies in a bin. ^dThe sSFR assigned to each stack is calculated using the median SFR and median ${M}_{* }$ of galaxies in each bin.

Download table as:  ASCIITypeset image

The ${M}_{* }$ stacks are shown in the SFR–${M}_{* }$ and sSFR–${M}_{* }$ planes in Figure 1. We determine the mean z ∼ 2.3 sSFR–${M}_{* }$ relation by fitting a linear function to the ${M}_{* }$ stacks in the sSFR–${M}_{* }$ plane. The best-fit coefficients are listed in Table 2. We determine the residuals about the mean ${M}_{* }$–sSFR relation by subtracting this mean relation from the sSFR estimate for each galaxy in the z ∼ 2.3 sample. In Figure 1, we compare our mean relation to that obtained by Shivaei et al. (2015) using the first 2 years of MOSDEF data. The two relations are generally consistent, although we find a slightly higher normalization, likely due in part to different SED fitting assumptions and the resulting stellar masses. We note that the MOSDEF z ∼ 2.3 Hβ detection limit biases the SFR in the lowest-mass bin high, and thus the fit we derive here may be steeper in the sSFR–${M}_{* }$ plane and shallower in the SFR–${M}_{* }$ plane than the true mean relation of the z ∼ 2.3 star-forming population. However, our results do not change significantly if we instead fit mean relations excluding the lowest-mass bin. Incompleteness in the lowest-mass bin due to the Hβ detection limit thus does not affect our conclusions.

Table 2.  Best-fit Linear Coefficients to z ∼ 2.3 Galaxy Properties as a Function of Stellar Mass

Property	Slope	Intercept
Star Formationa
SFR	0.67	−5.33
sSFR	−0.33	−5.33
Line Ratiosa
O3N2	−0.94	10.72
N2	0.59	−6.82
N2O2	0.41	−5.03
O32	−0.56	5.51
O3	−0.34	3.90
R23	−0.12	2.03
12+log(O/H)b
O3N2	0.30	5.30
N2	0.34	5.01
N2O2	0.30	5.27
O32	0.26	5.79

Notes.

^aCoefficients for $\mathrm{log}(X)=m\times \mathrm{log}({M}_{* }/{M}_{\odot })+b$, where X is the appropriate galaxy property, m is the slope, and b is the intercept. ^bCoefficients for $12+\mathrm{log}({\rm{O}}/{\rm{H}})=m\times \mathrm{log}({M}_{* }/{M}_{\odot })+b$, where m is the slope, b is the intercept, and 12+log(O/H) is determined using the corresponding line ratio.

Download table as:  ASCIITypeset image

In Figure 2, we present the excitation- and metallicity-sensitive line ratios O3N2, N2, N2O2, and O32 versus ${M}_{* }$. The line ratios O3N2 and O32 are sensitive to the ionization parameter, containing both a high and low ionization energy ionic species. N2 is also sensitive to changes in ionization parameter, as well as the nitrogen abundance (N/H), such that higher ionization parameter and lower N/H lead to lower N2. In contrast, N2O2 is primarily sensitive to changes in the N/O abundance ratio (Kewley & Dopita 2002). We observe a clear progression toward higher ionization parameter and lower N/O and N/H at fixed stellar mass from z ∼ 0 to z ∼ 2.3. These results are in agreement with what has been found for other samples at z > 2 (e.g., Erb et al. 2006; Steidel et al. 2014; Shapley et al. 2015; Holden et al. 2016). While recent results for high-redshift samples have suggested the possibility of redshift evolution of physical conditions such as the ionizing spectrum, ionization parameter, or N/O at fixed O/H (e.g., Kewley et al. 2013a; Masters et al. 2014, 2016; Steidel et al. 2014, 2016; Shapley et al. 2015; Sanders et al. 2016b; Strom et al. 2017), the observation of both higher ionization parameter and lower nitrogen abundance at fixed ${M}_{* }$ is difficult to explain without chemical evolution (i.e., lower O/H at fixed ${M}_{* }$ with increasing redshift) playing a major role. We determine the mean relation between these strong-line ratios and ${M}_{* }$ by fitting linear functions to the ${M}_{* }$ stacks, and present the best-fit coefficients in Table 2.

Zoom In Zoom Out Reset image size

All four line ratios indirectly trace oxygen abundance, since ionization parameter is tightly anticorrelated with metallicity (Dopita et al. 2006a, 2006b; Pérez-Montero 2014; Sánchez et al. 2015), while N/O and N/H are correlated with metallicity due to the secondary production channel of nitrogen (Kewley & Dopita 2002; Pettini & Pagel 2004; Pérez-Montero & Contini 2009). The observation of such correlations in the local universe has led to the construction of metallicity calibrations based on these strong-line ratios. The translation of the strong-line ratios O3N2, N2, N2O2, and O32 into metallicity is shown by the right set of y-axes in Figure 2. All four panels are matched to the same range in metallicity (7.7 < 12 + log(O/H) < 8.9). Under the assumption that metallicity calibrations do not strongly evolve with redshift, the presence of higher excitation and lower N/O at z ∼ 2.3 corresponds to lower metallicity at fixed ${M}_{* }$ compared to z ∼ 0. We find that z ∼ 2.3 galaxies are offset toward lower metallicities at fixed ${M}_{* }$ by 0.37, 0.25, 0.46, and 0.25 dex on average for metallicities based on O3N2, N2, N2O2, and O32, respectively. The best-fit linear MZR for these indicators can be found by passing the line-ratio versus ${M}_{* }$ fits through the metallicity calibrations (Equations (7), (8), (10), and (11)). We list the best-fit coefficients for the MZR in Table 2, and calculate O/H residuals relative to these best-fit linear relations.

The ability of such locally calibrated relations to accurately estimate nebular metallicities in high-redshift galaxies has been called into question by recent studies (e.g., Steidel et al. 2014, 2016; Sanders et al. 2015; Shapley et al. 2015; Strom et al. 2017). We note here that the method we employ to search for a ${M}_{* }$–SFR–Z relation at z ∼ 2 using residuals around mean relations is immune to changes in the normalizations of metallicity calibrations since we are looking at relative differences in metallicity, but would be affected by changes in the slopes. The interpretation of our results could also be affected by changes in ionized gas physical conditions that correlate with sSFR at fixed ${M}_{* }$. We discuss these potential systematic effects on our results in Section 4.1.

As previously mentioned, while it is of interest to determine whether or not z ∼ 0 and high-redshift galaxies lie on the same ${M}_{* }$–SFR–Z relation (i.e., the FMR holds out to high redshifts), it is more fundamental to establish whether the metallicities of high-redshift galaxies display a secondary dependence on SFR at fixed ${M}_{* }$. The existence of a high-redshift ${M}_{* }$–SFR–Z relation is a prerequisite for interpreting the position of high-redshift galaxies relative to the FMR observed at z ∼ 0. We look for a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 by searching for correlated scatter around the MZR and ${M}_{* }$–sSFR relation.

In the left column of Figure 3, we show the z ∼ 2.3 residuals around the mean ${M}_{* }$ relations for the line ratios O3N2, N2, and N2O2, and the corresponding O/H residuals, all as a function of residuals around the z ∼ 2.3 mean ${M}_{* }$–sSFR relation. We show both z ∼ 2.3 individual galaxies and the ${M}_{* }$–ΔsSFR stacks, the latter color-coded by stellar mass. There is a clear trend among both the stacks and individual galaxies such that the sSFR residuals correlate with O3N2 residuals and anticorrelate with residuals in N2 and N2O2. These trends are as expected in the presence of a ${M}_{* }$–SFR–Z relation, given that O3N2 decreases with increasing metallicity while N2 and N2O2 increase. Furthermore, the correlations are statistically significant based on Spearman correlation tests. The correlation coefficients and p-values are shown in Table 3. The corresponding trends in Δlog(O/H) are similar for all three metallicity indicators: Δlog(O/H) decreases with increasing Δlog(sSFR). This trend in the residuals demonstrates that z ∼ 2.3 galaxies with higher SFR have lower O/H at fixed stellar mass, confirming the existence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3. This is the first time that such a relation has been clearly demonstrated to exist at this redshift.

Zoom In Zoom Out Reset image size

We quantify the strength of the z ∼ 2.3 ${M}_{* }$–SFR–Z relation by fitting linear relations to the median Δlog(O/H) in bins of Δlog(sSFR). The best-fit lines are shown in the middle column of Figure 3, and the best-fit slopes are given in Table 3. These slopes are inconsistent with a flat relation (i.e., no dependence of O/H on SFR at fixed ${M}_{* }$) at 3–4σ. The ${M}_{* }$–ΔsSFR stacks fall on the same relations as the medians of the individual z ∼ 2.3 galaxies. The high-mass ${M}_{* }$–ΔsSFR stacks most clearly follow the same relation as the medians of the individual points, while the low-mass ${M}_{* }$–ΔsSFR stacks are noisier but all fall within 1σ of the best-fit relations. The position of the stacks relative to the individual detections suggests that galaxies with non-detections in these emission-line ratios lie on the same anticorrelations. We conclude that requiring line-ratio detections in the individual galaxies does not significantly bias these results.

The right column of Figure 3 displays the line-ratio and O/H residuals as a function of Δlog(sSFR) for the z ∼ 0 ${M}_{* }$–SFR stacks, color-coded by ${M}_{* }$. In each panel, we include the best-fit line to the z ∼ 2.3 sample for comparison. The z ∼ 0 stacks show similar trends at log(M_*/M_⊙) ≲ 10.0, but display the well-known decrease in the strength and disappearance of the SFR dependence at high stellar masses (Mannucci et al. 2010; Yates et al. 2012; Salim et al. 2014), as evidenced by the tendency of high-mass stacks to lie closer to Δlog(O/H) = 0 at fixed Δlog(sSFR/yr⁻¹). A weakening of the z ∼ 0 sSFR dependence at Δlog(sSFR/yr⁻¹) < 0 is also apparent in all three indicators, as pointed out by Salim et al. (2014). There is no evidence for such a flattening of the ${M}_{* }$–SFR–Z relation at high masses or Δlog(sSFR/yr⁻¹) < 0 among the z ∼ 2.3 ${M}_{* }$–ΔsSFR stacks. However, we note that the mass range in which the z ∼ 0 sample displays little to no SFR dependence (log(M_*/M_⊙) ≳ 10.5) lies above the mass range probed by the z ∼ 2.3 sample.

We show the line-ratio and MZR residual plots for O32 in Figure 4. There is no clear trend toward lower metallicity with increasing Δlog(sSFR) among the z ∼ 2.3 galaxies when utilizing the O32 indicator. There is no statistically significant correlation present (Table 3). In contrast, a clear ${M}_{* }$–SFR–Z relation is present at z ∼ 0 for O32. We discuss the lack of SFR dependence at z ∼ 2.3 for the O32 MZR in Section 4.3.

Zoom In Zoom Out Reset image size

The strong-line ratio O3 is sensitive to both ionization parameter and metallicity. As such, trends indicative of a ${M}_{* }$–SFR–Z relation should be present in the O3 residuals as well. We show O3 as a function of ${M}_{* }$ in the left panel of Figure 5. The best-fit coefficients are given in Table 2. The right panel of Figure 5 displays the O3 residuals as a function of sSFR residuals for individual z ∼ 2.3 galaxies and the ${M}_{* }$–ΔsSFR stacks. A clear correlation is present, with Spearman correlation coefficient r_s = 0.47 and probability of being drawn from an uncorrelated distribution of 1.6 × 10⁻¹³. The correlation between Δlog(O3) and Δlog(sSFR) is consistent with a decrease in metallicity (increase in O3) as sSFR increases at fixed ${M}_{* }$. This result agrees with what is found using O3N2, N2, and N2O2.

Zoom In Zoom Out Reset image size

The strong-line ratio R23 as a function of ${M}_{* }$ and the Δlog(R23) versus Δlog(sSFR) deviation plot are presented in Figure 6. R23 is primarily sensitive to O/H because it is a ratio of both low and high ionization state oxygen lines to a hydrogen Balmer line, but has significant ionization parameter dependence as well (Kewley & Dopita 2002). R23 is not a strong function of ${M}_{* }$ for the z ∼ 2.3 sample, which is anticipated since the R23 values lie close to the turnover regime in the R23–O/H relation, as described in Section 2.4. However, R23 does increase slightly with decreasing ${M}_{* }$, suggesting that our sample mostly lies on the upper metal-rich R23 branch where R23 increases with decreasing O/H. In this case, a positive correlation between Δlog(R23) versus Δlog(sSFR) is expected if a ${M}_{* }$–SFR–Z relation is present. This positive correlation is observed in the right panel of Figure 6 with high statistical significance (r_s = 0.60, p-value = 7.3 × 10⁻¹⁸), in agreement with the trends observed using O3N2, N2, N2O2, and O3.

Zoom In Zoom Out Reset image size

One important question is whether the SFR dependence we find in the z ∼ 2.3 MZR is a result of selection effects. We investigated the effect of requiring Hβ S/N ≥ 3 on these results by repeating the analysis to search for a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 without this requirement. Because we cannot reliably correct for reddening without detections of both Hβ and Hα, we utilized SFRs estimated from FAST SED fitting instead of dust-corrected Hα luminosities. Furthermore, this test could only be performed for the portions of the analysis that involve line ratios that do not require reddening correction (i.e., O3, N2, and O3N2). We found similar results based on these line ratios, although the trends were somewhat weaker for a few reasons. First, there is significant scatter in the correlation between SFR(SED) and SFR(Hα) (Reddy et al. 2015). Second, SFR(Hα) traces star formation on <10 Myr timescales while SFR(SED) probes longer timescales of ∼100 Myr. It is thus expected that signatures of a ${M}_{* }$–SFR–Z relation are weaker when relying on SFR(SED), since the ${M}_{* }$–SFR–Z relation is thought to be driven by a reaction of SFR and metallicity on short timescales to recent accretion of metal-poor gas (Mannucci et al. 2010; Davé et al. 2012). Finally, a smaller dynamic range of SFR at fixed ${M}_{* }$ is probed in this test because ${M}_{* }$ and SFR(SED) are not independent of one another, both being estimated by FAST SED fitting. This covariance artificially tightens the ${M}_{* }$–(s)SFR relation, making it more difficult to resolve SFR dependence of metallicity at fixed ${M}_{* }$. Because our results persist when utilizing SFR(SED) without the Hβ detection criterion, we conclude that requiring Hβ S/N ≥ 3 does not strongly bias our results.

Having established the existence of a ${M}_{* }$–SFR–Z relation for z ∼ 2.3 star-forming galaxies, we next address the question of whether or not these galaxies lie on the z ∼ 0 ${M}_{* }$–SFR–Z relation. In other words, we test whether or a FMR exists that extends out to z ∼ 2.3. Many studies have attempted to address this question by extrapolating parameterizations of the z ∼ 0 FMR (e.g., Lara-López et al. 2010; Mannucci et al. 2010) to the high-SFR and high-sSFR regime occupied by high-redshift galaxies (e.g., Christensen et al. 2012; Wuyts et al. 2012, 2014; Belli et al. 2013; Henry et al. 2013; Cullen et al. 2014; Yabe et al. 2014; Zahid et al. 2014b; Kashino et al. 2017). Such extrapolations may not be representative of the true z ∼ 0 ${M}_{* }$–SFR–Z relation because a parametric form must be assumed that may not represent the underlying physical relation, and the small fraction of z ∼ 0 objects with SFRs similar to those of high-redshift galaxies will not carry much weight toward the fit, allowing the possibility of a poor fit in the high-SFR regime. Recent studies have shown the benefit of performing non-parametric comparisons instead, relying on the small number of z ∼ 0 objects with extreme SFRs to directly compare at fixed ${M}_{* }$ and SFR (Salim et al. 2015; Sanders et al. 2015; Ly et al. 2016).

We directly compare the position of the z ∼ 2.3 ${M}_{* }$–ΔsSFR stacks and the z ∼ 0 ${M}_{* }$–SFR stacks in the mass–metallicity plane in Figure 7, with metallicities based on the O3N2, N2, and N2O2 indicators. Both high- and low-redshift stacks are color-coded by SFR on the same scale so that a direct comparison of metallicity at fixed ${M}_{* }$ and SFR is possible. The highest-SFR z ∼ 2.3 stack, with log(M_*/M_⊙) = 10.2 and log(SFR/M_⊙ yr⁻¹) ≈ 2.0, does not have any local analogue in the AM13 ${M}_{* }$–SFR sample. However, all other z ∼ 2.3 stacks have analogous z ∼ 0 counterparts matched in ${M}_{* }$ and SFR. In the O3N2, N2, and N2O2 ${M}_{* }$–Z planes, z ∼ 2.3 galaxies display metallicities that are systematically lower than their z ∼ 0 counterparts by ∼0.1 dex at fixed ${M}_{* }$ and SFR. This offset is consistent across a range of ${M}_{* }$ and Δlog(sSFR) using three different metallicity indicators. We found a similar offset using only O3N2 and N2 with a smaller z ∼ 2.3 MOSDEF sample in Sanders et al. (2015). We conclude that there is not a FMR that can simultaneously match the properties of star-forming galaxies from z ∼ 0 out to z ∼ 2.3.

Zoom In Zoom Out Reset image size

We note that the presence of a z ∼ 2.3 ${M}_{* }$–SFR–Z relation can be seen in Figure 7 in the ${M}_{* }$–ΔsSFR stacks, with the lowest-SFR bins at a given ${M}_{* }$ falling above the mean z ∼ 2.3 relation and the highest-SFR bins falling below. This trend is most evident in the O3N2 and N2O2 mass–metallicity planes (Figure 7) for the z ∼ 2.3 stacks with log(M_*/M_⊙) = 10.0 − 10.3, as shown by a progression from higher to lower O/H as SFR increases, and the symbol color changes from blue to purple to pink. The presence of such SFR gradients perpendicular to the mean z ∼ 2.3 MZR is a manifestation of the sSFR dependence shown in Figure 3.

We directly compare the strong-line ratios O3, O32, and R23 at fixed ${M}_{* }$ and SFR for z ∼ 0 and z ∼ 2.3 stacks in Figure 8. The z ∼ 2.3 stacks display higher O3 and R23 than z ∼ 0 stacks at fixed ${M}_{* }$ and SFR, suggesting higher excitation and lower metallicity at fixed ${M}_{* }$ and SFR at z ∼ 2.3. This result is consistent with what was found using O3N2, N2, and N2O2. The z ∼ 2.3 ${M}_{* }$–ΔsSFR stacks also display separation perpendicular to the mean O3–${M}_{* }$ and R23–${M}_{* }$ relations as a function of SFR, consistent with the existence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3. The O32 values of the z ∼ 2.3 stacks are not obviously offset from the z ∼ 0 stacks matched in ${M}_{* }$ and SFR. Once again, results based upon O32 are inconsistent with those from O3N2, N2, N2O2, O3, and R23. We discuss the differences in results based on O32 in Section 4.3.

Zoom In Zoom Out Reset image size

We have clearly demonstrated the existence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 based on multiple emission-line ratios (O3N2, N2, N2O2, O3, and R23) and metallicities inferred using z ∼ 0 calibrations. We furthermore showed that z ∼ 2.3 galaxies are offset ∼0.1 dex lower in metallicity from z ∼ 0 galaxies at fixed ${M}_{* }$ and SFR, arguing against the existence of a redshift invariant ${M}_{* }$–SFR–Z relation. We now consider whether the appearance of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 and its evolution with respect to z ∼ 0 could be produced by other effects in the absence of metallicity variations. Recent work has suggested that the physical conditions of ionized gas in star-forming regions evolve with redshift (e.g., Steidel et al. 2014, 2016; Shapley et al. 2015; Sanders et al. 2016b; Kashino et al. 2017; Strom et al. 2017). The physical conditions of interest are the ionization parameter, N/O abundance ratio, electron density, and the shape of the ionizing spectrum. Changes in these properties lead to changes in emission-line ratios at fixed nebular oxygen abundance (Kewley et al. 2013a). The ionization parameter, N/O ratio, and shape of the ionizing spectrum must be considered at fixed metallicity because of correlations between these properties and the nebular and stellar metallicity.

The only evolving physical condition for which there is a consensus is the electron density. The typical electron density in star-forming regions increases with redshift, and is an order of magnitude higher at z ∼ 2 compared to z ∼ 0 (e.g., Steidel et al. 2014; Sanders et al. 2016b; Kashino et al. 2017). The typical densities at z ∼ 0 and z ∼ 2.3 (25 and 250 cm⁻³, respectively) are both well below the density where strong-line ratios become strongly affected (∼1000 cm⁻³; Steidel et al. 2014). Thus, evolution in electron density has a negligible impact on our results.

There is not a consensus about the evolution (or lack thereof) of N/O, ionization parameter, and hardness of the ionizing spectrum at fixed metallicity. While this ambiguity makes it difficult to assess the impact of such evolution on metallicity studies, the MOSDEF sample crucially provides access to multiple strong-line ratios, enabling us to determine if evolution in a particular property leads to false conclusions. It has been suggested that high-redshift galaxies have elevated N/O at fixed O/H compared to z ∼ 0 galaxies, explaining observed evolution in line-ratio diagrams involving nitrogen (Masters et al. 2014; Steidel et al. 2014; Shapley et al. 2015; Sanders et al. 2016b). The elevated N/O could be caused by an unusually high occurrence rate of Wolf–Rayet stars (Masters et al. 2014) or a dilution of O/H by pristine gas inflows while N/O remains roughly constant (Köppen & Hensler 2005; Amorín et al. 2010; Sanders et al. 2016b). Other studies have argued for a harder ionizing spectrum at fixed nebular abundance in high-redshift star-forming regions (Steidel et al. 2016; Strom et al. 2017). Some have also suggested a higher ionization parameter at fixed metallicity as sSFR increases, leading to an elevated ionization parameter in high-redshift galaxies due to their larger sSFR than z ∼ 0 galaxies on average (Kewley et al. 2013a, 2013b, 2015, 2016; Kashino et al. 2017). In what follows, we consider the possible effects of these proposed scenarios on our results.

4.1.1. Nitrogen-to-Oxygen Ratio

4.1.2. Ionization Parameter and Hardness of the Ionizing Spectrum

The existence of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 confirms that the current theoretical framework for galaxy growth through the interplay of inflows, outflows, and star formation is applicable at high redshifts. A ${M}_{* }$–SFR–Z relation is predicted to exist at high redshifts in analytical chemical evolution models (Finlator & Davé 2008; Davé et al. 2012; Lilly et al. 2013), as well as cosmological hydrodynamical simulations including chemical evolution (Ma et al. 2016; Davé et al. 2017; De Rossi et al. 2017). There are both observational evidence and theoretical predictions that the ${M}_{* }$–SFR–Z relation is a manifestation of a more fundamental relation with gas mass, where the SFR is modulated by the amount of gas in a galaxy (Bothwell et al. 2013, 2016a, 2016b; Hughes et al. 2013; Zahid et al. 2014a; Ma et al. 2016; Davé et al. 2017). Confirmation of a ${M}_{* }$–SFR–Z relation at z ∼ 2.3 suggests that galaxy chemical evolution is linked to gas content at high redshifts as well. The increasing number of observations of the cold gas content of high-redshift galaxies (e.g., Tacconi et al. 2013) will allow for investigations of the relation between ${M}_{* }$, SFR, Z, and gas fraction at high redshifts.

By comparing residuals in metallicity and sSFR around the MZR and mean ${M}_{* }$–sSFR relation, we have quantified the strength of the SFR dependence of the ${M}_{* }$–SFR–Z relation at both z ∼ 0 and z ∼ 2.3. Current cosmological hydrodynamical simulations predict that the strength of the SFR dependence does not change with redshift (Davé et al. 2017; De Rossi et al. 2017). Using the mufasa simulations (Davé et al. 2016), Davé et al. (2017) predicted the slope of the Δlog(O/H) versus Δlog(sSFR) relation, finding the preferred slope to be −0.16 and independent of redshift. We have quantified this measure of the ${M}_{* }$–SFR–Z relation for our z ∼ 2.3 sample, providing best-fit slopes based on the O3N2, N2, and N2O2 metallicity indicators in Table 3.

We compare our observed z ∼ 2.3 Δlog(O/H) versus Δlog(sSFR) relations to the prediction of Davé et al. (2017) in the middle column of Figure 3, and perform the same comparison for the z ∼ 0 sample in the right column. In each panel, the solid black line is the best-fit z ∼ 2.3 relation and the dashed black line shows the prediction of Davé et al. (2017). We find weaker SFR dependence at z ∼ 2.3 than predicted by Davé et al. (2017) based on O3N2 and N2, and stronger SFR dependence when using N2O2 to estimate metallicites. The O3N2, N2, and N2O2 best-fit slopes all agree with the Davé et al. (2017) prediction within 2σ. Taking the uncertainty-weighted average of the three z ∼ 2.3 slopes yields a slope of −0.14 ± 0.022, which is within 1σ of a slope of −0.16. Our results are thus consistent with the predictions of the mufasa simulations, although both the measurement and systematic uncertainties remain large.

The SFR dependence of the z ∼ 2.3 sample displays a similar strength to that of z ∼ 0 galaxies with log(M_*/M_⊙) <10.0 and Δlog(sSFR/yr⁻¹) > 0 when using the N2 and N2O2 indicators (Figure 3). When using the O3N2 indicator, the z ∼ 0 relation appears to be somewhat stronger than at z ∼ 2.3. Overall, the observations do not suggest large changes in the strength of the SFR dependence of metallicity at fixed ${M}_{* }$ and SFR over the past 10.5 Gyr. We caution that a detailed comparison of the strength of the z ∼ 0 and z ∼ 2.3 ${M}_{* }$–SFR–Z relation requires a more complete understanding of the redshift evolution of gas physical conditions and the corresponding effects on metallicity indicators. Recalibration of strong-line ratio metallicity relations at high redshift using electron-temperature-based direct metallicities provides a promising avenue to eliminate these systematics (e.g., Jones et al. 2015; Sanders et al. 2016a).

The ∼0.1 dex offset toward lower metallicity at fixed ${M}_{* }$ and SFR observed from z ∼ 0 to z ∼ 2.3 demonstrates that the ${M}_{* }$–SFR–Z relation is not redshift invariant. While this shift in metallicity is fairly small, a systematic offset is observed across more than an order of magnitude in ${M}_{* }$ and SFR using four different metallicity-sensitive line ratios (O3N2, N2, N2O2, and O3). Observing the evolution of the ${M}_{* }$–SFR–Z relation in multiple line ratios simultaneously allowed us to show that this metallicity offset cannot be a false inference due to potential evolution of metallicity calibrations with redshift.

In Sanders et al. (2015), we found a similar evolution in O/H at fixed ${M}_{* }$ and SFR, and speculated that such an offset could be introduced if the inflow rate exceeded the sum of the SFR and outflow rate such that the gas reservoir grows faster than it can be processed via star formation. While such a gas-accumulation phase is predicted to exist, both models and observations suggest that it should only occur at z > 4 (Papovich et al. 2011; Davé et al. 2012). This phase can be made to reach z ∼ 2 in models with no outflows (Davé et al. 2012), which are clearly unphysical given the observational constraints on the occurrence of galactic winds at high redshifts (e.g., Steidel et al. 2010). It is thus unlikely that we are observing the buildup of galaxy gas reservoirs.

If the most fundamental relation is between ${M}_{* }$, metallicity, and gas content and that relation is redshift invariant, then an evolving ${M}_{* }$–SFR–Z relation would be indicative of an evolving relation between SFR and gas mass. If a non-evolving relation between ${M}_{* }$, gas fraction, and metallicity exists, then galaxies at z ∼ 0 and z ∼ 2.3 have the same gas mass at fixed ${M}_{* }$ and metallicity. The observed trend of decreasing metallicity at fixed ${M}_{* }$ and SFR with increasing redshift could then be understood as decreasing SFR at fixed ${M}_{* }$ and gas fraction (i.e., lower star formation efficiency). However, an inferred decrease in z ∼ 2 star formation efficiency at fixed ${M}_{* }$ is in conflict with the interpretation of recent observations of cold gas in high-redshift galaxies. Scoville et al. (2016, 2017) found that at fixed ${M}_{* }$, star formation efficiency per unit ISM gas mass increases with increasing redshift using ALMA observations of luminous star-forming galaxies at z ∼ 1−6. Furthermore, molecular gas depletion timescales of galaxies on the mean ${M}_{* }$–SFR relation are shorter at z ∼ 1−3 than at z ∼ 0 (Tacconi et al. 2013; Genzel et al. 2015), suggesting an increase in star formation efficiency with redshift. Reconciling these observations with the observed evolution of the ${M}_{* }$–SFR–Z relation requires an evolving relationship among ${M}_{* }$, gas fraction, and metallicity such that metallicity decreases at fixed ${M}_{* }$ and gas fraction with increasing redshift. Since a decreasing star formation efficiency with increasing redshift at fixed ${M}_{* }$ conflicts with observational constraints, we must investigate other drivers of ${M}_{* }$–SFR–Z evolution.

The cause of an evolving ${M}_{* }$–SFR–Z relation can also be understood through analytical chemical evolution models. The gas-regulator model of Lilly et al. (2013) allows for an evolving ${M}_{* }$–SFR–Z relation if the gas consumption timescale and/or the mass-loading factor (η = outflow rate/SFR) at fixed ${M}_{* }$ change with redshift. In this context, an evolving relation toward lower metallicity at fixed ${M}_{* }$ and SFR with increasing redshift suggests that, at fixed ${M}_{* }$, either the gas consumption timescale is longer (i.e., star formation efficiency is lower) at higher redshifts or η is higher at earlier times. The equilibrium model of Davé et al. (2012) similarly predicts lower O/H at fixed ${M}_{* }$ and SFR at high redshifts if η increases at fixed ${M}_{* }$ with increasing redshift. Although a lower star formation efficiency at higher redshift appears to be in conflict with the observations detailed above, a larger typical high-redshift mass-loading factor is at least broadly consistent with observational constraints over a wide range of ${M}_{* }$ and redshift (e.g., Pettini et al. 2002; Rupke et al. 2005; Weiner et al. 2009; Steidel et al. 2010; Heckman et al. 2015; Chisholm et al. 2017). We do note, however, that estimates of η carry large systematic uncertainties.

In both the equilibrium model of Davé et al. (2012) and the gas-regulator model of Lilly et al. (2013), the ISM metallicity at fixed ${M}_{* }$ and SFR decreases if the metallicity of the inflowing gas decreases. Although accreted gas is relatively unenriched compared to the ISM, it is not metal-free. Numerical simulations have shown that accreted gas is a mix of relatively pristine gas from the IGM and enriched gas that was ejected in past outflows (e.g., Oppenheimer et al. 2010; Segers et al. 2016; Anglés-Alcázar et al. 2017). Recycled gas can even become the dominant source of fuel for star formation at z < 1 (Oppenheimer et al. 2010). Thus evolution toward lower metallicity at fixed ${M}_{* }$ and SFR at z ∼ 2.3 could also be explained by a decrease in the metallicity of infalling gas at high redshift caused by a lower importance of recycled gas accretion relative to pristine accretion from the IGM. Numerical simulations agree with this scenario. Indeed, in the models of Davé et al. (2011), an increase in the metallicity of infalling gas relative to that of the ISM was the primary driver of the increase in metallicity at fixed ${M}_{* }$ with decreasing redshift. Anglés-Alcázar et al. (2017) found that the relative importance of recycled gas to total accretion increases with decreasing redshift in the FIRE zoom-in simulations and additionally suggest that transfer of metals between galaxies via galactic winds occurs at z = 0, further increasing the metallicity of infalling gas at z ∼ 0 compared to that at z ∼ 2. We conclude that the offset toward lower O/H at fixed ${M}_{* }$ and SFR with increasing redshift is most likely caused by an increase in the mass-loading factor at fixed ${M}_{* }$ from z ∼ 0 to z ∼ 2 and a decrease in the metallicity of infalling gas due to a higher relative importance of pristine accretion from the IGM at z ∼ 2.

While results based on O3N2, N2, N2O2, O3, and R23 are all consistent with respect to both the existence of a z ∼ 2.3 ${M}_{* }$–SFR–Z relation and the ∼0.1 dex offset toward lower metallicity at fixed ${M}_{* }$ and SFR, performing the same analysis using O32 yields different results. The Δlog(O32) versus Δlog(sSFR) diagram shows no evidence for SFR dependence of O/H at fixed ${M}_{* }$ (Figure 4). Furthermore, there is little if any evolution in O32 at fixed ${M}_{* }$ and SFR from z ∼ 0 to z ∼ 2.3 (Figure 8, middle panel). The disagreement between results based on O32 compared to those based on the other line ratios requires closer examination.

Since the emission lines in O32 are widely separated in wavelength, it is natural to suspect that incorrect reddening corrections may give rise to discrepant results based on O32. If there are large galaxy-to-galaxy deviations from the assumed Cardelli et al. (1989) attenuation curve, a trend in O32 with SFR at fixed ${M}_{* }$ could be washed out. However, systematic effects associated with galaxy-to-galaxy variations in the attenuation curve do not appear to have a significant effect on the O32 results, since we see strong trends in N2O2 and R23 (which also require a reddening correction) with SFR at fixed ${M}_{* }$.

SFR dependence of O32 at fixed ${M}_{* }$ could also be hidden if our assumed dust-correction recipe overestimates E(B–V)_gas for galaxies lying above the mean z ∼ 2.3 SFR–${M}_{* }$ relation and/or underestimates E(B–V)_gas in the opposite regime. In this case, O32 values at low Δlog(sSFR) would be overestimated while the opposite would occur at high Δlog(sSFR), masking the presence of a ${M}_{* }$–SFR–Z relation. Interestingly, such a scenario would steepen Δlog(R23) versus Δlog(sSFR) and Δlog(N2O2) versus Δlog(sSFR) beyond what is expected based on the ${M}_{* }$–SFR–Z relation because R23 increases and N2O2 decreases with increasing SFR at fixed ${M}_{* }$ (Figures 3 and 6). Indeed, we find that the slope of Δlog(O/H) versus Δlog(sSFR) is steepest when metallicities are based on N2O2, and that Δlog(R23) versus Δlog(sSFR) displays the strongest correlation of all line ratios considered here. It is therefore possible that our dust-correction recipe leads to either underestimated E(B–V)_gas for z ∼ 2.3 galaxies below the SFR–${M}_{* }$ relation, overestimated E(B–V)_gas for objects above the SFR–${M}_{* }$ relation, or both. If true, a bias in E(B–V)_gas dependent on SFR suggests that the nebular attenuation curve systematically changes with SFR at fixed ${M}_{* }$ at z ∼ 2.3. Studies of the stellar continuum attenuation curve and the differential reddening between the continuum and nebular emission lines at high redshifts have similarly suggested a systematic dependence of the shape of the attenuation curve on SFR and sSFR (Kriek & Conroy 2013; Price et al. 2014; Reddy et al. 2015). We note that the resulting bias would affect the measured strength of the SFR dependence of O/H at fixed ${M}_{* }$ using O3N2, N2, and O3 as well, since biases in E(B–V)_gas will change Δlog(sSFR).

Biases in reddening corrections may also explain why no evolution in O32 at fixed ${M}_{* }$ and SFR is observed from z ∼ 0 to z ∼ 2.3. If E(B–V)_gas is systematically overestimated in the z ∼ 2.3 sample, then the high-redshift O32 values would be underestimated and fall closer to the z ∼ 0 galaxies, reducing any offset present. If this is the case, then the offset in O/H at fixed ${M}_{* }$ and SFR based on N2O2 is overestimated. In Figure 7, the metallicity offset at fixed ${M}_{* }$ and SFR based on N2O2 is larger than what is found using O3N2 and N2, consistent with an overestimate of the reddening correction at z ∼ 2.3. It is thus plausible that biases in the dust-correction can explain the discrepant O32 results.

Unlike O3, O3N2, N2, N2O2, and R23, O32 directly probes the ionization parameter and is only sensitive to metallicity because of the anticorrelation between ionization parameter and metallicity. This anticorrelation arises because gas-phase metallicity and stellar metallicity are correlated, giving rise to harder ionizing spectra and, consequently, higher ionization parameters in lower metallicity star-forming regions (Dopita et al. 2006a, 2006b). The O32 results could be explained if ionization parameter does not depend on metallicity at z ∼ 2.3. However, such a scenario would weaken trends in O3N2 and O3, both of which are strongly sensitive to ionization parameter in addition to metallicity. We observe clear trends in O3N2 and O3 as a function of SFR at fixed ${M}_{* }$. Furthermore, O32 is clearly still related to metallicity given that O32 is anticorrelated with stellar mass, and a similar MZR is obtained when using O32 compared to those based on O3N2, N2, and N2O2 (Figure 2). We conclude that ionization parameter and metallicity are still anticorrelated at z ∼ 2.3, though not necessarily in the same form as at z ∼ 0.

Another possibility is that the scatter in the ionization parameter versus metallicity relation is significantly larger at high redshifts than at z ∼ 0. If this is the case, additional scatter would be introduced in the relationship between O32 and metallicity, potentially washing out trends in O32 with SFR at fixed ${M}_{* }$. This scenario can be tested by measuring the scatter in O32 at fixed metallicity at z ∼ 2.3, where metallicity is based on temperature-sensitive auroral-line measurements, but a data set on which to perform such an analysis does not yet exist.

Based on the information currently available, biases arising from the reddening correction are the most likely explanation of the discrepant O32 results. These biases are such that E(B–V)_gas is systematically overestimated at z ∼ 2.3, and the magnitude of this overestimate increases with increasing SFR at fixed ${M}_{* }$.