A CENSUS OF OXYGEN IN STAR-FORMING GALAXIES: AN EMPIRICAL MODEL LINKING METALLICITIES, STAR FORMATION RATES, AND OUTFLOWS

We draw our local sample from the SDSS DR7, which consists of ∼900,000 galaxies spanning a redshift range of 0 < z < 0.7 (Abazajian et al. 2009). The survey has a Petrosian limiting magnitude of r_P = 17.8 covering 8200 deg². The spectra have a nominal spectral range of 3900–9100 Å and a spectral resolution of R ∼ 2000. We make use of the ugriz-band photometry available for each object (Stoughton et al. 2002) and the publicly released emission-line fluxes measured by the MPA-JHU group.⁵ We refer to the sample presented here as the "metallicity-selected SDSS sample."

We correct for dust extinction in the emission lines by inferring a reddening correction from the Balmer decrement. For case B recombination with electron temperature T_e = 10⁴ K and electron density n_e = 10² cm⁻³, the intrinsic Hα/Hβ ratio is expected to be 2.86 (Osterbrock 1989). We get the intrinsic color excess, E(B − V), and the correction for dust attenuation using the extinction law of Cardelli et al. (1989) and a corresponding R_v = 3.1. We note that the results of this study are not dependent on our choice of a particular extinction law.

We select a pure star-forming sample of local emission-line galaxies from the SDSS DR7. We first distinguish star-forming galaxies from active galactic nuclei (AGNs) by constraining the ionizing radiation source using the [O iii] λ5007, [N ii] λ6584, Hβ, and Hα emission lines (Baldwin et al. 1981; Veilleux & Osterbrock 1987; Kauffmann et al. 2003; Kewley et al. 2006). In particular, we remove galaxies using the equation given in Kewley et al. (2006) where

$$\begin{equation} \mathrm{log([{\rm O\,\mathsc{iii}}]/H\beta)} > 0.61/(\mathrm{log([{\rm N\,\mathsc{ii}}]/H\alpha)} - 0.05 ) + 1.3. \end{equation} \tag{ 1 }$$

In order to avoid aperture effects, we require a g-band fiber aperture covering fraction >30% in addition to imposing a lower redshift limit of 0.04 (Kewley et al. 2004). The median covering fraction for the metallicity-selected SDSS sample is 38%. Kewley et al. (2006) find that the SDSS sample is incomplete at higher redshifts, and in order to minimize evolutionary effects, we also impose an upper limit redshift cutoff of z = 0.1.

In order to establish comparable samples, galaxies in the local sample are selected from SDSS using the same selection criteria as the DEEP2 sample. In particular, galaxies are selected to have signal-to-noise ratio (S/N) of Hβ > 3, σ_R23 < 2, and equivalent width of Hβ > 4 Å. Here, σ_R23 is the error in the R23 parameter, which is the ratio of the oxygen nebular emission ([O ii] λ3727 doublet and [O iii] λλ4959, 5007) to Hβ. These particular selection criteria give us a sample of ∼20,000 star-forming galaxies in the limited redshift range of 0.04 < z < 0.1.

Our sample of intermediate-redshift star-forming galaxies is taken from the DEEP2 survey (Davis et al. 2003). Details of sample selection and properties are given in Zahid et al. (2011); here we summarize the data selection. The survey consists of ∼45,000 galaxies targeted mostly in the redshift range of 0.7 > z > 1.4 by applying a color preselection using BRI-band photometry (Coil et al. 2004). The survey has a limiting magnitude of R_AB = 24.1 and covers 3.5 deg². The spectra have a nominal spectral range of 6500–9100 Å and a resolution of R ∼ 5000. The emission-line equivalent widths are measured in Zahid et al. (2011), and we adopt the same values here.

The sample selection is based on the spectral and photometric properties of the galaxies. We reduce AGN contamination by first requiring that log(R23) < 1 (L. J. Kewley et al. 2012, in preparation). Using the color separation for blue and red galaxies parameterized by Willmer et al. (2006), Weiner et al. (2007) conclude that only a small fraction of blue galaxies in DEEP2 appear to harbor AGNs, whereas a large fraction of red emission-line galaxies show evidence of AGN emission. We further limit AGN contamination by removing 48 galaxies in the sample that are classified as red galaxies using the color division of Willmer et al. (2006).

Given the nominal spectral coverage, the redshift range of galaxies in our sample is limited by the necessity to simultaneously observe both the [O ii] λ3727 doublet and the [O iii] λ5007 emission lines, which are required for chemical analysis. We further require that the S/N Hβ > 3, the error of the R23 emission-line ratio, σ_R23, be less than 2, and the equivalent width of Hβ > 4 Å. Finally, due to ambiguity in the metallicity determination of DEEP2 galaxies at low stellar masses (see Section 3.2), we also remove all galaxies with M_* < 10^9.2 M_☉. This selection criterion gives us a sample of 1348 star-forming galaxies in the limited redshift range of 0.75 < z < 0.82.

Erb et al. (2006a, E06 hereafter) determine the MZ relation from 87 star-forming galaxies at z ∼ 2.2 selected from a larger sample of 114 galaxies described in Erb et al. (2006b). The galaxies are selected on the basis of their rest-frame UV colors, and redshifts are determined from their UV spectra. All galaxies in the sample have UGRK-band photometry. Most also have J-band photometry, and 32 galaxies have been observed at 3.6, 4.5, 5.4, and 8.0 μm with IRAC on board the Spitzer Space Telescope. Hα spectra were obtained using NIRSPEC on the Keck II telescope.

The metallicity for these galaxies is determined from the emission-line ratio of [N ii] λ6584 to Hα. The [N ii] λ6584 line is sensitive to metallicity, with the strength of the line decreasing with decreasing metallicity. The S/N of the individual galaxy spectra is insufficient to measure the weak [N ii] λ6584 line. In order to increase the S/N of their spectra and to increase the chance of detecting [N ii] emission line at low metallicities, E06 stack 14 or 15 individual galaxy spectra binned by stellar mass into 6 composite spectra. The [N ii] λ6584 and Hα emission-line fluxes and metallicities are measured from these composite spectra.

We use the Le Phare⁶ code developed by S. Arnouts and O. Ilbert to estimate the galactic stellar mass. This code estimates the stellar masses of galaxies by comparing photometry with stellar population synthesis models in order to determine the mass-to-light ratio, which is then used to scale the observed luminosity (Bell et al. 2003b; Fontana et al. 2004). We synthesize magnitudes from the stellar templates of Bruzual & Charlot (2003) and use a Chabrier (2003) initial mass function (IMF). The 27 models have two metallicities and seven exponentially decreasing star formation models (SFR ∝e^−t/τ) with τ = 0.1, 0.3, 1, 2, 3, 5, 10, 15, and 30 Gyr. We apply the extinction law of Calzetti et al. (2000) allowing E(B − V) to range from 0 to 0.6, and the stellar population ages range from 0 to 13 Gyr. The median statistical error for the derived stellar masses, determined from propagating the uncertainty in the photometry, is 0.15 dex. Though systematic effects may be large (Drory et al. 2004; Conroy et al. 2009), we have consistently measured the stellar masses for our different samples, giving us a robust relative measure (Swindle et al. 2011; R. Swindle 2012, private communication).

We account for emission-line contributions by taking the Kennicutt (1998a) relation between the synthesized UV luminosity and SFR and the emission lines. This treatment accounts for Hα, Hβ, [O ii] λ3727, and [O iii] λλ4959, 5007 (Ilbert et al. 2009). In our sample, making no correction for the emission-line contributions does not significantly alter our mass determinations. We adopt the median of the mass distribution and take the 68% confidence interval as a measure of the error. In Zahid et al. (2011), we compare this method with the method used by the MPA/JHU group to determine stellar masses of the SDSS galaxies. We find that the two estimates differ by a constant offset of ∼0.2 dex and that the dispersion between the two methods is 0.14 dex.

E06 measure stellar masses using a similar method of comparing photometry with stellar population synthesis models (Erb et al. 2006c). However, an important difference is that they measure the "total" stellar mass, which is the integral of the SFR over the lifetime of the galaxy. They find that this stellar mass is ∼10%–40% higher than the instantaneous (or what they term "current living") stellar mass. However, most mass estimates found in the literature are generally the instantaneous and not "total" stellar mass. We recalculate the stellar masses for their individual galaxies and bin the data according to their original binning. In Figure 1, we compare the mass estimates for each of their six bins. The x-axis shows our mass estimate, and the y-axis is the original mass estimate from E06. The error bar represents the rms dispersion of stellar masses in each mass bin. For the five higher mass bins, we find that there is a constant offset such that the E06 estimates are 0.29 dex higher than our estimate. For the lowest mass bin, there appears to be a near one-to-one agreement, though the 0.29 dex offset is within the errors. This is most likely due to the fact that the lowest mass galaxies are younger and therefore have a much closer agreement between their instantaneous and "total" stellar mass. Moreover, we do not concern ourselves too much with the lowest mass bin as the metallicity measure is only an upper limit for this bin. Since we do not rebin the E06 data in determining the metallicity, we adopt their stellar mass values but subtract 0.29 dex to make them consistent with our estimates.

Zoom In Zoom Out Reset image size

Figure 2(a) shows the distribution of the stellar masses for all three samples.

Zoom In Zoom Out Reset image size

We use two strong-line methods to determine the metallicity of galaxies in our sample. In this section, we only present the parameterization of the calibration and defer a detailed discussion of the methods until Section 8.3. For the SDSS and DEEP2 sample, we determine metallicities using the calibration of Kobulnicky & Kewley (2004). This method relies on the R23 and O32 parameters, which are defined as

$$\begin{equation} R23 = \frac{\mathrm{[{\rm O\,\mathsc{ii}}]}\,\lambda 3727 + \mathrm{[{\rm O\,\mathsc{iii}}]} \,\lambda\lambda 4959, 5007}{\mathrm{H}\beta } \end{equation} \tag{ 2 }$$

and

$$\begin{equation} O32 = \frac{\mathrm{[{\rm O\,\mathsc{iii}}]}\,\lambda\lambda 4959, 5007}{\mathrm{[{\rm O\,\mathsc{ii}}]}\,\lambda 3727}. \end{equation} \tag{ 3 }$$

Here, the ratio implies the ratio of the measured line intensities. We have used the assumption that the ratio of the fluxes of [O iii]λ5007 to [O iii]λ4959 is 3 (Osterbrock 1989). Due to the higher S/N of the [O iii]λ5007 line, for the sum of the [O iii]λ4959 and [O iii]λ5007 flux we adopt a value of 1.33 times the [O iii]λ5007 flux.

The R23 strong-line method is known to be sensitive to the ionization parameter. The ionization parameter characterizes the ionization state of the gas and quantitatively represents the number of ionizing photons per second per unit area divided by the hydrogen density. The ionization parameter can be constrained by measuring the ratio of the line intensity of the same element at two ionization stages. The method of Kobulnicky & Kewley (2004) does this explicitly by using the O32 line ratio. Because both the metallicity and ionization parameter are interdependent, an iterative scheme is used, the details of which are provided in the Appendix of Kewley & Ellison (2008). In particular, we use Equations (A4) and (A6).

The DEEP2 data are not flux calibrated, so we use the line equivalent widths with a correction applied for Balmer absorption (Zahid et al. 2011). For consistency, we determine the SDSS metallicities using line equivalent widths as well. Several studies have established that line equivalent widths can be substituted for line fluxes when measuring metallicities if data are not flux calibrated or a reliable reddening estimate is unavailable (Kobulnicky & Phillips 2003; Moustakas et al. 2010; Zahid et al. 2011). In particular, we test this on our SDSS sample and find that the dispersion between the metallicities measured using equivalent widths and dereddened line fluxes is ∼0.05 dex, which is less than intrinsic uncertainties of the strong-line method. The average difference in the MZ relation derived using the two methods is <0.01 dex.

One issue with using the R23 method is that metallicity is not a monotonic function of R23. For a particular value of the R23 parameter there are two possible values of metallicity known as the upper and the lower branch metallicity. This degeneracy is generally alleviated by employing a second metallicity diagnostic that relies on line ratios that are monotonic with metallicity. For our SDSS sample of galaxies, the use of [N ii]/Hα reveals that the vast majority (∼99%) of galaxies are on the upper branch. For our DEEP2 sample, the nominal wavelength coverage of the spectra does not allow us to simultaneously observe the R23 lines and [N ii]/Hα. We assume that all galaxies with M_* > 10^9.2 M_☉ lie on the upper metallicity branch (for more details see Figure 5 in Zahid et al. 2011).

The calibration of Kobulnicky & Kewley (2004) is based on theoretical photoionization models. Empirical methods rely on calibrating strong-line ratios to metallicities determined using temperature-sensitive auroral lines. E06 determine metallicities using the empirical calibration of Pettini & Pagel (2004). The metallicity is given by

$$\begin{equation} \mathrm{12+log(O/H)} = 8.90 + 0.57 \times N2, \end{equation} \tag{ 4 }$$

where N2 = log([N ii]λ6584/Hα). Using ∼28, 000 galaxies from the SDSS DR4, Kewley & Ellison (2008) derive constants for converting between various diagnostics. We use these conversion constants to consistently determine metallicities for all three samples using both the theoretical and empirical calibrations of Kobulnicky & Kewley (2004) and Pettini & Pagel (2004).

The metallicity of a galaxy is traditionally defined as a number density of oxygen relative to hydrogen and is given as 12 + log(O/H). For this study, we require the mass density of oxygen relative to hydrogen. We convert number density to mass density using the relation

$$\begin{eqnarray} \mathrm{log(}Z_g) &= & 12 + \mathrm{log(O/H)} - 12 - \mathrm{log} \left(\frac{M_O/M_H}{XM_H + YM_{He}} \right) \nonumber \\ & = & \mathrm{log(O/H)} - \mathrm{log} \left(\frac{15.999/1.0079}{0.75\times 1.0079 + 0.25\times 4.0026} \right)\!. \nonumber\\ \end{eqnarray} \tag{ 5 }$$

For the remainder of the paper, Z_g refers to the gas-phase mass density of oxygen relative to hydrogen. Also, when referring to metallicity or gas-phase oxygen abundance, we mean to refer to Z_g.

Figure 2(b) shows the distribution of the metallicities for all three samples.

The SFRs for the SDSS sample are derived by the MPA/JHU group using the technique of Brinchmann et al. (2004) with additional improvements given by Salim et al. (2007). The strong emission lines of each galaxy are fit using the nebular emission models of Charlot & Longhetti (2001). All emission lines contribute to constraining the dust attenuation, though the largest contribution comes from the Hα/Hβ ratio. At the median redshift of the SDSS, about 1/3 of the galaxy light is contained within the 3'' fiber. In order to account for losses, Brinchmann et al. (2004) derive an aperture correction, which was improved upon by Salim et al. (2007). We convert from a Kroupa to Chabrier IMF by dividing by 1.06.

We estimate the SFRs for the DEEP2 and E06 sample from the Hβ and Hα emission lines, respectively. In order to reliably estimate the SFR, it is necessary to make a correction to the Balmer line flux due to dust extinction. Since the nominal wavelength coverage of both the DEEP2 and E06 spectra does not allow us to observe multiple order Balmer lines, we are not able to derive a correction from the Balmer decrement. Garn & Best (2010) have shown that for a sample of star-forming galaxies from the SDSS the extinction determined from the Hα/Hβ ratio can be predicted from the stellar mass, Hα luminosity, or metallicity of the galaxy with a dispersion of ∼0.1 dex. The relation between stellar mass and dust extinction found by Garn & Best (2010) is shown to be valid out to z ∼ 1.5 (Sobral et al. 2012), and a relation between dust extinction, stellar mass, and metallicity is also observed in galaxies at z ∼ 2 (Reddy et al. 2010).

Xiao et al. (2012) obtain a slightly better fit (dispersion of 0.07 dex) by incorporating galaxy inclination and using a different parameterization. We determine E(B − V) for our sample of SDSS galaxies from the Balmer decrement assuming the Cardelli et al. (1989) extinction law and parameterize E(B − V) as a function of stellar mass and metallicity using a similar formulation to Xiao et al. (2012). The color excess is given by

$$\begin{equation} E(B-V) = (p_{0} + p_{1} Z^{p_2}) \times M^{p_3}, \end{equation} \tag{ 6 }$$

where Z = 10^{(12 + log(O/H) − 8)} and M = M_*/10¹⁰. The metallicities are derived using the calibration of Kobulnicky & Kewley (2004). We perform a nonlinear least-squares fit using the MPFIT⁷ set of routines (Markwardt 2009). The data are weighted by their 1σ statistical uncertainties, and the errors are propagated to the parameters. The fit returns p₀ = 0.12 ± 0.01, p₁ = 0.041  ±  0.006, p₂ = 0.77  ±  0.06, and p₃ = 0.240  ±  0.002. The left panel of Figure 3 shows the mean E(B − V) binned by stellar mass and color-coded by metallicity. The fit to these binned data is shown by the color-coded curves. The right panel shows the residuals between the observed and fitted color excess. There is a slight tail to the distribution owing to the lower quality of the fit to the high stellar mass, high-metallicity, and high color excess data. The rms dispersion between the observations and the fit is ∼0.11 dex.

Zoom In Zoom Out Reset image size

We measure the SFRs for the DEEP2 sample of galaxies using the Hβ luminosity. In order to obtain a flux calibration, we integrate the DEEP2 spectra over the I-band filter response and compare to the observed I-band magnitude. This gives us an estimate for the conversion between counts on the detector and flux given in ergs s⁻¹. This procedure accounts for both throughput and slit losses assuming that the line-to-continuum emission ratio is constant inside and outside of the slit. This assumption is reasonable as slit losses are generally small because the vast majority of DEEP2 galaxies have effective radius smaller than the 1'' slit width (Weiner et al. 2007). We linearly fit the flux-calibrated continuum in a 80 Å window about the Hβ emission line and scale the equivalent width by the continuum to obtain an Hβ flux.

We use the extinction law of Cardelli et al. (1989) and the fitted E(B − V) values given by Equation (6) to correct the observed Hβ flux in the DEEP2 sample for dust attenuation. We obtain the luminosity from the flux using the relation L(Hβ) = F_Hβ × 4πD²_L, where F_Hβ is the dereddened Hβ flux and D_L is the luminosity distance determined from the redshift. We determine the SFR from the Kennicutt (1998a) relation between SFR and Hα luminosity given by

$$\begin{equation} \mathrm{SFR} = 7.9 \times 10^{-42} L(\mathrm{H}\alpha)\, [\mathrm{erg} \,\, \mathrm{s}^{-1}]. \end{equation} \tag{ 7 }$$

We take L(Hα) = 2.86 × L(Hβ) and scale the SFR from Salpeter to Chabrier IMF by dividing by 1.7.

E06 infer their SFRs from the Hα luminosity corrected for dust extinction using the E(B − V) determined from spectral energy distribution fitting. A factor of two aperture correction was applied, and the SFRs were determined using the Kennicutt (1998a) relation converted to a Chabrier IMF. We redetermine the SFRs for their 87 galaxies by applying an extinction correction determined from Equation (6) to their observed Hα luminosities. Following E06, we apply a factor of two aperture correction to the observed luminosity. We determine the SFR from Equation (7) and bin data according to E06.

Figure 2(c) shows the distribution of the SFRs for all three samples.

In Figure 4 we plot the MS relation for the SDSS (blue squares), DEEP2 (black circles), and E06 (red triangles) star-forming galaxies at redshifts of z = 0.07, 0.78, and 2.26, respectively. For the SDSS, DEEP2, and E06 samples the median SFR is determined in bins of stellar masses. For the SDSS, DEEP2, and E06 samples there are ∼700, ∼90, and ∼15 galaxies in each data bin. The flattening in the relation observed at lower stellar masses is due to incompleteness (see below) such that the lower SFR galaxies fall out of our selected samples at lower stellar masses. Down to log(M_*) ∼ 9.2 for the SDSS and DEEP2 sample and log(M_*) ∼ 9.5 for the E06 sample the relation between stellar mass and SFR is well described by a power law.

Zoom In Zoom Out Reset image size

For the SDSS and DEEP2 samples, the error bars plotted on the data points indicate the interval containing 68% of the data within each mass bin. For the E06 sample the error bars indicate the standard deviation of the distribution within each mass bin. The median scatter in the SFRs is 0.25, 0.27, and 0.24 dex for the SDSS, DEEP2, and E06 samples, respectively. The scatter in the three samples is comparable and is roughly constant with stellar mass. The observed scatter in our samples is comparable to the ≲ 0.3 dex found in studies examining the M_*–SFR relation out to z ∼ 2 (Salim et al. 2007; Elbaz et al. 2007; Noeske et al. 2007b; Daddi et al. 2007; Pannella et al. 2009).

We fit a linear relation between log(Ψ) and log(M_*/M_☉) in the mass range where the samples are monotonically declining using the linfit.pro routine in IDL. The form of the fit is given by

$$\begin{equation} \mathrm{log(SFR)} = \delta + \gamma \cdot [ \mathrm{log}(M_\ast /M_\odot) - 10], \end{equation} \tag{ 8 }$$

where δ is the logarithm of the SFR at 10¹⁰ M_☉ and γ is the power-law index. The dashed blue, solid black, and dotted red lines in Figure 4 are the fits to the SDSS, DEEP2, and E06 samples, respectively. We give the fit parameters and the errors in Table 2. The errors for the SDSS and DEEP2 fit parameters are bootstrapped. The errors for the E06 sample are determined from propagating the standard error of each data bin, $\sigma /\sqrt{n}$, where σ is the standard deviation of the distribution and n is the number of objects in each bin. The power-law index of the relation is consistent (within the small errors) for the SDSS and DEEP2 samples. We have applied a factor of two aperture corrections to all galaxies in the E06 sample. This correction is likely overcorrecting the lower mass galaxies, which could explain the shallower slope observed in the MS relation. The zero point of the MS relation shows evolution to higher values with increasing redshift.

In all three of our samples the MS relation flattens out at lower stellar masses (see Figure 4). The galaxies in all three of the samples are selected such that chemical properties can be determined from the spectra. Determining chemical properties from the emission lines of star-forming galaxies generally requires high- S/N observations, and therefore samples with determined metallicities are more susceptible to incompleteness. We demonstrate that the observed flattening is due to incompleteness by comparing our sample of ∼20,000 galaxies for which we have well-determined metallicities with a larger, more complete sample of SDSS star-forming galaxies.

We select the larger comparison sample by first removing AGNs using the [O iii]/Hβ versus [N ii]/Hα diagram (Kauffmann et al. 2003; Kewley et al. 2006) and imposing the similar redshift restrictions (z < 0.1) as our metallicity-selected SDSS sample. This gives us a parent sample of ∼200,000 star-forming galaxies. The aperture-corrected SFRs given in the SDSS DR7 rely on a Bayesian technique principally fitting Hα and Hβ lines to determine extinction-corrected SFRs (Brinchmann et al. 2004; Salim et al. 2007). We select galaxies from the parent sample by requiring an S/N of 5 in the Hα and Hβ emission lines. This gives us a sample of ∼140,000 star-forming galaxies with z < 0.1 and detections of Hα and Hβ. We note that the distribution of SFRs in our comparison sample is largely independent of choice of S/N cuts and therefore the observed incompleteness in our metallicity-selected SDSS sample likely results from our limit on the equivalent width of Hβ.

Figure 5 shows the SFRs for both our metallicity-selected SDSS sample of galaxies and the larger comparison sample. The black data points with error bars and dashed line are the same as in Figure 4. The solid black line is the median SFRs of our comparison sample sorted into 75 bins of stellar mass, and the gray curves are the 68% and 95% contours of the data. The MS relation in the larger comparison sample has a slope of 0.65, which is slightly more shallow than our metallicity-selected SDSS sample. This slope is consistent with the slope determined by Salim et al. (2007) for the MS relation of SDSS galaxies. The zero point of the larger sample is slightly lower as compared to our metallicity-selected SDSS sample. Over the range where we fit the MS relation, the difference in the SFRs for the two samples is <0.1 dex, and both relations show similar scatter.

Zoom In Zoom Out Reset image size

From this comparison, it is clear that the flattening observed in the MS relation in Figure 4 is due to incompleteness. The slightly higher zero point of our metallicity-selected SDSS sample suggests that there is a slight bias against low star-forming galaxies over all mass ranges. This is consistent with our expectation that the higher S/N and equivalent width of Hβ required in the spectra for chemical analysis select against galaxies with weak lines and low levels of star formation. Despite this selection effect, Figure 5 demonstrates that incompleteness does not significantly bias our determination of the MS relation.

Noeske et al. (2007b) examine the MS relation for DEEP2 galaxies in the redshift range of 0.2 < z < 1. They determine SFRs from the 24 μm flux and add in the SFRs determined from the DEEP2 emission lines with no extinction-correction to account for unobscured star formation. They determine the fit to the MS relation in the range where their sample is >95% complete. The MS relation determined by Noeske et al. (2007b; but see Dutton et al. 2010) at z = 0.78 is consistent with our determination for our DEEP2 sample at the same redshift, despite the fact that our sample is selected differently and we have determined SFRs from extinction-corrected emission lines. Our correction for dust extinction relies on the correlation between extinction, stellar mass, and metallicity (Section 4.3). This necessarily introduces a correlation between stellar mass and SFR. The fact that the MS determined from our DEEP2 sample is consistent within the errors to the relation fit by Noeske et al. (2007b) implies that our dust correction is not introducing a spurious correlation between stellar mass and SFR. Furthermore, we conclude that the MS relation we determine at z = 0.78 from our DEEP2 sample is not biased due to incompleteness.

The completeness of the E06 sample at z = 2.26 is much more uncertain. This sample is UV selected, which may bias against star-forming galaxies with low levels of star formation if they are present at this redshift. However, we note that the observed relation determined from the E06 sample is consistent within the errors to the zero point and scatter found by Daddi et al. (2007), using a highly complete sample at z ∼ 2, and Pannella et al. (2009), who use radio stacking to determine SFRs. This suggests that the observed MS relation in the E06 sample is not significantly biased.

With a well-defined MS relation at several redshifts, we can parameterize the SFR as a function of mass and redshift. We seek a parameterization where the variables of mass and redshift are separated such that the SFR is given by

$$\begin{equation} \Psi (M_\ast , z) = \psi _{o}(z) \cdot f(M_\ast). \end{equation} \tag{ 9 }$$

Here Ψ(M_*, z) is the SFR as a function of stellar mass and redshift, ψ_o(z) is the zero point, which evolves with redshift, and f(M_*) is the relation between stellar mass and SFR. This definition explicitly assumes that the slope of the relation between stellar mass and SFR has no time dependency.

As discussed in Section 4.1, the SFR as a function of stellar mass is well fit by a power law in stellar mass with an index near unity. We therefore parameterize the stellar mass dependency as

$$\begin{equation} f(M_\ast) = (M_\ast /10^{10} \,M_\odot)^{\gamma }, \end{equation} \tag{ 10 }$$

where the only free parameter, γ, is the power-law index of the relation. With this definition for f(M_*), the zero point, ψ_o(z), is the SFR of a 10¹⁰ M_☉ galaxy. The zero point of the relation as a function of redshift, as we show below, is well fit by an exponential, and we parameterize it as

$$\begin{equation} \psi _o(z) = \alpha \cdot \mathrm{exp}(\beta z). \end{equation} \tag{ 11 }$$

Here, α is the zero point at z = 0 and β determines the rate of exponential increase with redshift.

Figure 6 shows a compilation taken from the literature for the values of γ (bottom panel) and ψ_o (top panel) as defined in Equations (10) and (11), respectively. For comparison, data from this study are plotted as black points on these plots. We have converted the stellar masses and SFRs for the data shown in Figure 6 such that they are consistent with the Chabrier (2003) IMF. The data selection and method of SFR determination vary. In the local universe, Salim et al. (2007) derive SFRs by fitting the UV/optical broadband SED. Elbaz et al. (2007) adopt the Brinchmann et al. (2004) aperture-corrected SFRs given in the DR4. At 0.2 < z < 1, the measurements of Noeske et al. (2007b)⁸ have SFRs determined from a combination of uncorrected emission-line fluxes and 24 μm fluxes. At z ∼ 1, Elbaz et al. (2007) determine SFRs from the infrared luminosity inferred from the 24 μm flux. At z ∼ 2, Daddi et al. (2007) determine SFRs from dust-corrected UV luminosities and Pannella et al. (2009) determine SFRs from stacking analysis of the radio continuum.

Zoom In Zoom Out Reset image size

Despite the differences in methodology, our determination of the zero point, ψ_o (top panel Figure 6), is remarkably consistent with values found in the literature. In the local universe, we find a marginally higher zero point, which is due to incompleteness (see Section 4.2). At z = 0.78 we determine the same relation (within the errors) as Noeske et al. (2007b). At z = 2.26 the relation determined from the E06 sample is consistent with the higher redshift data. Out to z ∼ 1, the slope of the MS relation, γ (bottom panel of Figure 6), is consistent in most of the studies. At z > 1 the slopes show some variation, though the variation is only about ±35%. We note that in this study, our results are not sensitive to the adopted value of γ because the redshift evolution of the SFR is significantly greater than the magnitude of the variation introduced by a ±35% change in γ.

In order to derive Ψ(M_*, z), we need to determine values for the three free parameters, α, β, and γ. For the SDSS and DEEP2 samples, the slopes are consistent (within the errors), and we adopt a constant value of 0.7 for γ in Equation (10) (see also Whitaker et al. 2012). This is shown by the dashed line in the bottom panel of Figure 6. To determine α and β, we fit a line of the form

$$\begin{equation} \mathrm{log}[\psi _o(z)] = \mathrm{log}(\alpha) + 0.434 \, \beta z \end{equation} \tag{ 12 }$$

to our data. The fit is shown by the dashed line in the top panel of Figure 6. We note that we have fit only to our data, but the fit is more or less consistent with all the studies compiled from the literature. At z > 2.26 we assume that ψ_o is flat, noting that our results are not strongly dependent on this extrapolation because in our model most of the stellar mass growth in star-forming galaxies in the local universe occurs at z < 2.26 (see also Leitner 2012).

With both ψ_o(z) and f(M_*) determined, we can parameterize the SFRs of star-forming galaxies as a function of stellar mass and redshift by entering the fitted results back into Equation (9). The SFR as a function of stellar mass and redshift is given by

$$\begin{equation} \Psi (M_\ast , z) = 2.00 \cdot \mathrm{exp}(1.33 z) \left(\frac{M_\ast }{10^{10}}\right)^{0.7} [M_\odot \, \mathrm{yr}^{-1}]. \end{equation} \tag{ 13 }$$

In Figure 7, we show the fitted relation plotted over the observations. Though the fitted slope of the MS relation at z = 2.26 was shallower (see Table 2), the data are still reasonably well described by the relation given in Equation (13).

Zoom In Zoom Out Reset image size

The total amount of oxygen produced by the stars within a galaxy is simply given by

$$\begin{equation} M_T^o = P_o \, \frac{M_\ast }{1-R}. \end{equation} \tag{ 21 }$$

The quantity M_*/(1 − R) represents the total amount of gas converted into stars. Several simplifying assumptions go into Equation (21): (1) the IMF is constant, (2) gas is instantaneously returned, and (3) the oxygen yield is independent of metallicity. In Section 8, we discuss systematics and uncertainties associated with each of these assumptions.

The total mass of oxygen in the gas phase in local star-forming galaxies is given by M^o_g = Z_g M_g, where M^o_g is the mass of oxygen in the gas-phase. Unfortunately, direct gas mass measurements are available for only a few galaxies in our local sample. We therefore estimate the gas mass of local galaxies using the observed relation between gas fraction and stellar mass for star-forming galaxies. Peeples & Shankar (2011) compile gas mass estimates (H i + H₂, with correction for helium) from several sources in the literature and give the binned cold gas fraction as a function of stellar mass (see Table 2 in their paper). The values are compiled from the data sets of McGaugh (2005), Leroy et al. (2008), and Garcia-Appadoo et al. (2009).

In Figure 9 we plot the log of the gas fraction, f_g, as a function of stellar mass. The gas fraction is defined such that f_g = M_g/M_*.⁹ The stellar mass estimates of Peeples & Shankar (2011) are taken from the MPA/JHU catalog, and we make a 0.2 dex correction for consistency (Zahid et al. 2011). The error bars are the 1σ uncertainties on the mean of log(f_g). The dashed line is a fit to the relation given by

$$\begin{eqnarray} \mathrm{log}(f_g) &=& (-0.28 \pm 0.03)\nonumber\\ && - (0.55 \pm 0.03)\cdot [\mathrm{log}(M_\ast /M_\odot) - 10]. \end{eqnarray} \tag{ 22 }$$

The errors in the parameters come from propagating the uncertainties in the mean of log(f_g). When calculating gas mass, we add 0.1 dex to the fitted relation at all stellar masses to account for ionized gas (see Table 1.1 in Tielens 2005). Using this relation, the gas mass is given by M_g = f_gM_*. We note that the gas fractions given by Equation (22) are greater than values inferred from inversion of the Schmidt–Kennicutt relation (see Figure 2 in Peeples & Shankar 2011). Using gas fractions determined from inverting the Schmidt–Kennicutt relation would therefore lower the estimate of the mass of oxygen in the gas phase and increase the oxygen deficit presented in Section 8.4.

Zoom In Zoom Out Reset image size

We can straightforwardly estimate the mass of oxygen in the gas phase of local star-forming galaxies from observed quantities by combining the relation for f_g given in Equation (22) with the MZ relation for local star-forming galaxies shown in Figure 8. This is given by M^o_g = Z_gf_gM_*.

The most difficult component of our oxygen census to estimate is the oxygen locked up in stars. The total amount of oxygen produced by stars and the amount of oxygen found in the gas phase of galaxies can be inferred from the present-day physical quantities of the stellar mass and gas-phase abundance. However, the oxygen locked up in stars is an accumulated property that is dependent on the star formation and chemical history. The amount of oxygen locked up in stars between some initial redshift t_i and some later redshift t is analytically given by

$$\begin{equation} M^{o}_{\ast }(t>t_i) = (1 - R)\, \int ^t_{t_i} Z_g(t^\prime) \, \Psi (t^\prime) \,dt^\prime. \end{equation} \tag{ 23 }$$

We have parameterized both the SFR (see Section 4) and gas-phase oxygen abundance (see Section 7) as a function of stellar mass and time and are therefore able to empirically constrain the total amount of oxygen locked up in stars by simply carrying out a numerical integration of Equation (23).

We simultaneously determine the star formation and stellar mass history of a galaxy using Equation (19). We set M_{*, i} = 10⁶ M_☉ and carry the integral from t_i out to t = 12.5 Gyr.¹⁰ Higher values of t_i correspond to younger, lower stellar mass star-forming galaxies in the local universe. We show the stellar mass history of three galaxies (t_i = 2.4, 4.0, and 6.4 Gyr) in Figure 10(a) and star formation history in Figure 10(b) assuming R = 0.35. The stellar mass history as a function of redshift allows us to determine the metallicity as a function of redshift for each galaxy using the procedure described in Section 7. We plot the chemical history in Figure 10(c). Entering the star formation history, Ψ(t) (Figure 10(b)), and the chemical history, Z_g(t) (Figure 10(c)), into Equation (23), we can determine the mass of oxygen locked up in stars as a function of redshift. This is shown in Figure 10(d).

Zoom In Zoom Out Reset image size

Gallazzi et al. (2005) and Panter et al. (2008) examine the stellar-mass−stellar-metallicity relation determined from stellar population modeling of SDSS spectra. Both groups derive a consistent relation, which implies that the imprint of the MZ relation from earlier epochs is observed in stellar populations of galaxies today. In Figure 11, we compare the metallicity of stars from our model analysis to that observed by Panter et al. (2008) using ∼300, 000 galaxies in SDSS. We estimate the mass-weighted stellar metallicity from our model of the star formation and chemical histories of galaxies (see Figures 10(b) and (c), respectively). In particular, we compute

$$\begin{equation} Z_\ast (M_\ast) = \frac{\int _{t_i}^{t} \Psi (t^\prime) Z(t^\prime) dt^\prime }{\int _{t_i}^{t} \Psi (t^\prime)dt^\prime } = \frac{\int _{t_i}^{t} \Psi (t^\prime) Z(t^\prime) dt^\prime }{M_\ast }, \end{equation} \tag{ 24 }$$

where Z_* is the mass-weighted stellar metallicity as a function of current stellar mass, M_*. We compare directly to the mass-weighted stellar metallicity relation observed and parameterized by Panter et al. (2008; see Equation 1 and Figure 6). In Figure 11 the Panter et al. (2008) fitted relation is plotted by the solid black curve and the 16% and 84% contours of the distribution are shown by the dotted black curve.

Zoom In Zoom Out Reset image size

As can be seen from Equation (24), our estimate of the metal content of stars is dependent on the choice of gas-phase abundance diagnostic. In Figure 11, we plot the mass-weighted stellar metallicity for both the Kobulnicky & Kewley (2004, red dashed curve) and the Pettini & Pagel (2004, green dot-dashed curve) calibration. For reference the MZ relation from Tremonti et al. (2004) is plotted by the yellow curve, and for consistency we adopt the solar metallicity value of 12+log(O/H)_☉ = 8.87 used by Panter et al. (2008). Estimates from the two different abundance diagnostics converted into stellar metallicities using Equation (24) envelope the relation from Panter et al. (2008). In this study, we have adopted a metallicity calibration that is an average of the Kobulnicky & Kewley (2004) and the Pettini & Pagel (2004) calibrations (see Section 7). Within the systematic uncertainties, our estimate of stellar metallicities using Equation (24) and the observed stellar metallicities from Panter et al. (2008) are consistent.

In Figures 12(a)–(c) we plot the mass of oxygen for our three census components. We plot the total mass of oxygen produced (blue curve), the mass of oxygen in the gas phase (black curve), and the mass of oxygen locked up in stars (red curve) as a function of stellar mass for star-forming galaxies in the local universe. The solid, dotted, and dashed curves are used to illustrate the effect varying the parameters has on different components of the census. The three free parameters in our model are the return fraction, R, the oxygen yield, P_o, and a constant offset in metallicity (which we refer to as ΔZ_g) since the zero point varies depending on the calibration. Figures 12(a)–(c) demonstrate the dependencies of the various components of our census on the free parameters of the model. In Figures 12(a)–(c) we vary R, P_o, and ΔZ_g, respectively, while keeping the other two parameters fixed. The three values adopted for the varying parameter and the adopted values of the constant parameters are shown in (a)–(c). In Figure 12(a) we see that variations in R mainly affect the total amount of oxygen produced with a small effect on the amount of oxygen locked up in stars. Figure 12(b) shows that varying P_o only affects the total amount of oxygen produced, and Figure 12(c) demonstrates that varying ΔZ_g affects both the stellar and gas-phase components of the census but has no effect on the total amount of oxygen produced.

Zoom In Zoom Out Reset image size

For a closed-box model, the total mass of oxygen produced equals the mass of oxygen found in the gas phase and stars. However, if metals are lost from the system, the equality will not hold. Thus, we have

$$\begin{equation} \Delta M^o = M_T^o - M_g^o - M_\ast ^o, \end{equation} \tag{ 25 }$$

where, as before, M^o_T is the total mass of oxygen produced and M^o_g and M^o_* are the mass found in the gas phase and stars, respectively. ΔM^o, which we refer to as the oxygen deficit, represents the total amount of oxygen produced that is unaccounted for by the gas-phase and stellar components. In Figure 12(d) we plot the oxygen deficit, which is the same for the three cases shown in Figure 12(a)–(c). The variation in the parameters is degenerate such that the oxygen deficit depends only on the variation in the quantity ΔZ_g(1 − R)/P_o, effectively reducing our three free parameters to one free parameter. We note that the stellar yield, y, is given by y = P_o/(1 − R). Therefore, our model is only sensitive to the changes in the zero point of the metallicity relative to the stellar yield. Higher (lower) values of ΔZ_g/y correspond to a smaller (larger) oxygen deficit. To first order, this is true for all models of chemical evolution relying on these quantities.

If we assume that all oxygen in galaxies originates from the stars within the galaxy, the oxygen deficit cannot be negative because a negative value would imply that there is more oxygen in the gas and stars than was produced by the galaxy. A physical model of the oxygen deficit requires a combination of the free parameters (ΔZ_g(1 − R)/P_o) that yields ΔM^o ⩾ 0. The dotted line in Figure 12(d) gives a negative oxygen deficit at M_* ≲ 10¹⁰ M_☉, implying that this is an unphysical model.

The functional form of the oxygen deficit with respect to stellar mass varies according to ΔZ_g(1 − R)/P_o. However, a generic feature of the oxygen deficit is that it monotonically increases with stellar mass and is substantial at high stellar masses regardless of our particular choice of ΔZ_g(1 − R)/P_o. In the case where the oxygen deficit is substantial at all stellar masses (solid curve in Figure 12), we have ΔM^o∝M^α_* with a near-unity value for the exponent (α = 1.14). When the oxygen deficit is large at all stellar masses (solid curve in Figure 12), the oxygen deficit does not strongly depend on the shape of the MZ relation or the shape of the relation between gas fraction and stellar mass. At more moderate values of ΔZ_g(1 − R)/P_o examined in our models (dashed curve in Figure 12), the relation between the oxygen deficit and stellar mass is more complicated. In Section 9.1, we argue that a model where the oxygen deficit is large at all stellar masses and is nearly proportional to the stellar mass is most consistent with independent observations of oxygen in the halos of star-forming galaxies.

In Figure 13, we plot the fractional oxygen deficit. In the top panel we plot the logarithm of δ_o, which we define as

$$\begin{equation} \delta _o = \frac{\Delta M^o}{M_\ast ^o + M_g^o}. \end{equation} \tag{ 26 }$$

δ_o is the oxygen deficit relative to the total oxygen content of the galaxy. In the bottom panel we plot the logarithm of ρ_o, which we define as

$$\begin{equation} \rho _o = \Delta M^o/M_g^o. \end{equation} \tag{ 27 }$$

δ_o is an important quantity reflecting key physical processes governing chemical evolution of galaxies but cannot be determined for individual galaxies as the fraction of oxygen locked-up stars cannot be measured. ρ_o does not account for oxygen locked up in stars but can be determined from observable quantities and is therefore an important parameter for comparison with observations.

Zoom In Zoom Out Reset image size

The fraction of gas returned to the ISM is mass dependent and therefore sensitive to the particular choice of an IMF, with theoretical values ranging between ∼20% and 50%. We have adopted this range in Figure 12(a). The physical mechanism regulating mass loss in stars varies by stellar mass. In low-mass stars, stellar winds during the red giant branch and asymptotic giant branch are the dominant mechanism. In intermediate-mass stars, winds during the asymptotic giant branch phase dominate, and for high-mass stars, both stellar winds and SNe contribute to stellar mass loss (Jungwiert et al. 2001, and references therein).

The return fraction is important in our models in that it sets the total amount of star formation required for a galaxy to have a particular current stellar mass. The total amount of gas converted into stars, M_gs, in a galaxy with a current stellar mass of M_* is given by M_gs = M_*/(1 − R). A galaxy with R = 0.5 will, to first order, form 40% more stars than a galaxy with R = 0.3 to achieve the same stellar mass. Because the amount of oxygen produced is directly proportional to the total amount of star formation, there will be a commensurate increase in total oxygen production.

We have also adopted the instantaneous recycling approximation (IRA) introduced by Tinsley (1980). The IRA assumes that the gas is returned to the ISM immediately following star formation. It is shown that the IRA works well for high-mass stars whose lifetimes are short compared to star formation timescales. Oxygen is a primary element, mainly produced in Type II SNe. In our models, the lifetimes of massive stars responsible for oxygen enrichment and gas mixing timescales are significantly shorter than the timescales of star formation, and therefore the IRA is reasonable in regard to oxygen production. However, the gas return from Type II SNe accounts for 10%–20% of the gas returned to the ISM, and a larger fraction is returned on longer timescales from intermediate- and low-mass stars (see Jungwiert et al. 2001).

Examining various IMFs, Leitner & Kravtsov (2011) show that the majority of stellar mass loss occurs in the first 1–2 Gyr (cf. their Figure 1). Noeske et al. (2007b, among others) argue that the MS relation implies that for most star-forming galaxies the SFR gradually declines. The gradual decline in SFRs is a generic feature of models assuming the continuity condition that galaxies build up stellar mass evolving along the observed MS relation (e.g., Peng et al. 2010; Leitner & Kravtsov 2011; Leitner 2012; this work). The IRA can be reasonably applied in situations where the SFR is not changing rapidly with redshift because, though the current generations of stars being formed will return the majority of their gas in 1–2 Gyr, the total gas return will be an integration over previous generations. If the SFR is not rapidly varying, the return fraction will reach a relatively constant value within a few Gyr.

In the Appendix, we apply a time-dependent treatment for gas return in determining the oxygen deficit. We no longer assume the IRA, so that oxygen is not instantaneously recycled, nor is the return fraction a constant. We show that our results are not sensitive to time dependencies in the return rate and adopting a constant and instantaneous return of gas is valid.

The production of oxygen in our models is defined by the oxygen yield, P_o. The oxygen yield represents the total amount of newly synthesized oxygen per unit stellar mass of star formation. Following Henry et al. (2000), we formally define the integrated oxygen yield as

$$\begin{equation} P_o = \int _{m_l}^{m_u} m p_o(m) \phi (m)dm. \end{equation} \tag{ 28 }$$

Here, m_l and m_u are the lower and upper limits of the mass range of stars formed, respectively, p_o(m) is the mass fraction of oxygen synthesized by a star of stellar mass m, and ϕ(m) is the normalized IMF. Uncertainties in P_o arise from uncertainties in the stellar yields, the IMF, and the mass range of integration. Here, we discuss only the stellar oxygen yields. For a more complete review of uncertainties in stellar yields of various elements see Romano et al. (2010).

The dependence of the oxygen yield on the IMF is clear from examination of Equation (28). The relative weight of stars in different mass ranges assuming different IMFs is given by Romano et al. (2005, Table 1). The relative weight of 8–40 M_☉ stars is a factor of 1.6 higher for the Chabrier IMF as compared to the Salpeter IMF, highlighting the uncertainties in the oxygen yields associated with the assumption of the IMF. These uncertainties are further exacerbated by the possibility of a varying IMF. Kroupa & Weidner (2003) argue that the integrated galactic initial mass function (IGIMF), the stellar mass distribution function for all stars within a galaxy, is an integration of the stellar IMF over the embedded cluster mass function (ECMF). There exists an empirical relation between the ECMF and the observed SFR, leading to a dependence of the IGIMF on the observed SFR. Peeples & Shankar (2011) show several models for the stellar yield, taking into consideration these variations.

Almost all oxygen in the universe is produced in massive stars (m ≳ 8 M_☉). Using detailed models of Type II SN explosions, Woosley & Weaver (1995) provide the nucleosynthetic yields for stars in the mass range of 11–40 M_☉ and for a range of metallicities. Thielemann et al. (1996) present similar model calculations, but in the mass range of 13–25 M_☉. Thomas et al. (1998) compare these two models, finding good agreement between the oxygen yields except at high stellar masses, where the yields calculated by Woosley & Weaver (1995) saturate. Henry et al. (2000) calculate the integrated oxygen yield using the models of Maeder (1992) and Woosley & Weaver (1995) and the Thielemann et al. (1996) models updated to include contributions from 40 and 70 M_☉ stars (Nomoto et al. 1997). Assuming a Salpeter IMF and a range of upper and lower masses, Henry et al. (2000) find that the oxygen yields range between 0.004 and 0.016 (see their Table 2). Using more sophisticated models, Kobayashi et al. (2006) calculate nucleosynthetic yields by assuming a range of metallicities, explosion energies (SNe and hypernovae), and metallicity-dependent mass loss. They calculate the integrated oxygen yield using a Salpeter IMF and a lower and an upper mass limit of 0.07 and 50 M_☉, respectively. They find that the oxygen yield is weakly dependent on metallicity and ranges between 0.006 and 0.008. The effect of oxygen yield variations on the estimate of the total amount of oxygen produced is degenerate with variations in the stellar mass determination. For example, Zahid et al. (2011) show that stellar masses estimated for the same sample of galaxies vary by −0.2 dex depending on methodology. In this study, we have taken a conservative approach and adopted the lower stellar mass estimate. Adopting the higher stellar mass estimate would have the same effect on the total amount of oxygen produced as increasing the oxygen yield by the same factor. In Figure 12(b), we adopt a range of 0.0056–0.0091 for the oxygen yield.

The uncertainty in the absolute nebular abundance calibration is a long-standing problem in observational astronomy. Kewley & Ellison (2008) show that the metallicity can vary by as much as 0.7 dex when using different abundance diagnostics for the same set of galaxies. This poses serious problems for understanding galactic chemical evolution. The uncertainties are largely due to the method of calibration used in establishing the diagnostic. Empirical methods rely on calibrating strong-line ratios using H ii regions and galaxies with metallicities determined using the so-called direct method, which uses temperature-sensitive auroral lines (e.g., Pettini & Pagel 2004). Empirical methods have several limitations, and theoretical calibrations of strong-line ratios using photoionization models have also been developed (e.g., Kewley & Dopita 2002; Kobulnicky & Kewley 2004; Tremonti et al. 2004). Detailed discussions of strengths and weaknesses of the various methods and calibrations can be found in Kewley & Ellison (2008) and Moustakas et al. (2010); we summarize the salient points here.

Peimbert (1967) points out that temperature fluctuations within nebular regions may lead to overestimates of the temperature. Temperature fluctuations are thought to be more problematic in metal-rich H ii regions, where efficient line cooling may lead to temperature inhomogeneities and strong temperature gradients (Garnett 1992). These overestimates of nebular temperatures lead to underestimates of the metallicity with the direct method (Stasińska 2002, 2005; Bresolin et al. 2006), and temperature inhomogeneities among the H ii regions within a galaxy also pose problems for metallicity determinations based on global spectra (Kobulnicky et al. 1999). A second and perhaps more pernicious problem is that the distribution of electrons within an H ii region may not follow a Boltzmann distribution but rather may be κ-distributed (a well-known particle distribution in plasma physics). A κ-distribution of electrons may account for the discrepant temperatures and metallicities inferred from auroral lines (Nicholls et al. 2012).

Empirical calibrations require observations of extremely faint auroral lines. In particular, the [O iii]λ4363 line is on the order of 100 times weaker than the strong optical oxygen lines (e.g., [O ii]λ3727, [O iii]λ5007). Furthermore, the line strength diminishes with increasing metallicity, due to efficient line cooling, and is only observed in low-metallicity H ii regions (≲ 0.5 Z_☉). The high-S/N spectra required to observe the auroral lines taken together with temperature inhomogeneities within individual H ii regions suggest that samples used to calibrate empirical methods may be biased and unreliable, particularly at higher metallicities.

Theoretical methods rely solely on photoionization models to calibrate strong-line ratios. These methods are not susceptible to observational limitations imposed by empirical calibrations, and metallicities are well constrained and parameterization is well defined over the full range of observed line ratios. The absolute calibration of metallicities is model dependent, and uncertainties are subject to the simplifying assumptions of the model. Indirect evidence suggests that these methods systematically overestimate the metallicity (e.g., Bresolin et al. 2009a; Kudritzki et al. 2012).

The metallicities in our samples are determined using either the theoretical diagnostic of Kobulnicky & Kewley (2004; SDSS and DEEP2) or the empirical diagnostic of Pettini & Pagel (2004; E06). These two methods differ by a constant offset of ∼0.3 dex. For our metallicity diagnostic, we therefore adopt the average of the metallicities determined using these two methods.

Kewley & Ellison (2008) have shown that the relative estimates of the metallicity are robust for various diagnostics. However, the absolute calibration varies significantly, and the shape of the MZ relation can change depending on the calibration adopted (see Figure 2 of Kewley & Ellison 2008). In our models, the oxygen deficit is sensitive to the shape of the MZ relation only in the case when the oxygen deficit is small. As we show in the following section, our model showing a large oxygen deficit (the solid line in Figure 12(d)) is most consistent with observations of oxygen in the halos of star-forming galaxies. Furthermore, the shapes of the MZ relations determined from the calibration of Kobulnicky & Kewley (2004) and Pettini & Pagel (2004) are consistent within the errors (see Figure 8) despite the fact that the two diagnostics are calibrated independently using a theoretical and empirical calibration, respectively.

Some uncertainty in the absolute metallicity scale is associated with depletion of metals onto dust. Several authors have inferred ∼0.1 dex depletion of oxygen onto dust grains (Jenkins 2004; Cartledge et al. 2006; Peimbert & Peimbert 2010). However, in NGC 300 the agreement between the gradient inferred from A and B supergiants and that from emission lines suggests that there is very little if any dust depletion in this galaxy (Bresolin et al. 2009a). It is likely that dust depletion ranges between 0 and 0.1 dex and perhaps is dependent on physical properties of galaxies such as metallicity or stellar mass.

The oxygen deficit we estimate is sensitive to the gas-phase abundance only at small oxygen deficits. If we assume a constant relative level of depletion (e.g., ∼0.1 dex), the effects of dust depletion on the oxygen deficit are more pronounced in lower mass galaxies, which have a smaller oxygen deficit. However, low-mass galaxies may also have lower levels of depletion as NGC 300¹¹ indicates. By comparing to independent observations of the mass of oxygen found in the halos of star-forming galaxies, we argue in Section 10.1 that a model with a large oxygen deficit is favored. In the case of a large oxygen deficit (solid curve in Figure 12), assuming that oxygen is depleted onto dust by 0.1 dex changes the oxygen deficit by 0.7 dex at a stellar mass of 10⁹ M_☉ and by 0.15 dex at a stellar mass of 10¹¹ M_☉. This level of dust depletion in our model with a large oxygen deficit does not alter the conclusions of this paper. Given the large uncertainties on how dust depleted oxygen is in the H ii region and systematic uncertainties of our metallicity calibration, we make no explicit correction for this effect. In Figure 12(c), we vary the adopted metallicity by [−0.2, +0.1] dex.

It is well established that star-forming galaxies in the local universe have spatial abundance gradients such that the central regions of galaxies are more metal-rich than their outskirts (e.g., Zaritsky et al. 1994). In our estimate for the amount of oxygen in the gas phase (see Section 8.2), we adopt an oxygen abundance determined from nebular emission within the 3 arcsec SDSS fiber aperture. We determine the oxygen abundance of local galaxies from data selected to have a covering fraction of >30%, which is shown to be sufficient for reproducing the global abundance (Kewley et al. 2005). We reproduce the global luminosity-weighted abundance, whereas an unbiased method for estimating the gas-phase metal content of galaxies would be to use a global metallicity estimator that reflects the average total gas-phase abundance (i.e., a gas-mass-weighted abundance). However, given current observations, it is not possible to assess what biases are introduced by adopting a luminosity-weighted abundance for estimating gas-phase metal content. Surveys of abundance gradients in nearby galaxies currently under way along with a new generation of radio instruments capable of measuring the gas content of large samples of galaxies will likely resolve this issue.

We speculate that the luminosity-weighted oxygen abundance is likely to overestimate the gas-mass-weighted abundance. In star-forming galaxies where luminosity is dominated by young stars, the luminosity-weighted abundance reflects the abundance of star-forming gas. The Schmidt–Kennicutt relation for star formation relates the star formation surface density to the gas surface density by Σ_SFR∝Σ^1.4_gas (Kennicutt 1998b). Because of the nonlinear relation between the two quantities, the luminosity-weighted metallicity will not be consistent with the gas-mass-weighted metallicity but instead will be biased toward higher gas surface density regions, which tend to lie in the central, metal-rich regions of galaxies. Therefore, the average metallicity of the gas is likely to be lower than the inferred luminosity-weighted abundance used in this study for estimating the oxygen mass in the gas phase. However, given that observations comparing the gas-mass- and luminosity-weighted abundances are currently not feasible, we do not attempt to make a correction for this effect.

Galactic outflows are a key physical process governing galaxy evolution (Veilleux et al. 2005). In the nearby universe they are commonly observed in starburst and post-starburst galaxies (e.g., Rupke et al. 2005; Martin 2006; Tremonti et al. 2007; Rich et al. 2010; Tripp et al. 2011) and are ubiquitous in higher redshift star-forming galaxies (e.g., Shapley et al. 2003; Weiner et al. 2009; Rubin et al. 2010; Steidel et al. 2010). The outflowing gas and metals may remain gravitationally bound to galaxies, in which case the material should be found to reside in the hot halos, or it may escape the galaxy potential well altogether.

The oxygen deficit may be partially resolved by oxygen present in the outer disks and hot halos of galaxies. Several studies have found flat abundance gradients in the outer disks of star-forming galaxies (Bresolin et al. 2009b; Werk et al. 2010, 2011; Bresolin et al. 2012). The amount of oxygen observed in the outer disks cannot be reconciled with the low levels of star formation, and transport of metals from the inner disk is the most likely scenario explaining the flat abundance gradients though the physical mechanism for the transport of metals to outer disks is not clearly understood. Oxygen has also been observed in absorption in the halos of star-forming galaxies (Chen & Mulchaey 2009; Prochaska et al. 2011). In a survey of 42 galaxies conducted using the Cosmic Origins Spectrograph (COS) on board the Hubble Space Telescope, Tumlinson et al. (2011) report ubiquitous detection of O vi in the halos of star-forming galaxies. They conclude that the mass of heavy elements and gas in the CGM may far exceed what is found within the galaxies themselves. From the column density of O vi, they estimate that the CGM of 10^9.5–10^11.5 M_☉ galaxies (out to 150 kpc) contains at least ∼10⁷ M_☉ of oxygen, where they have assumed an ionization correction factor of 0.2 to account for oxygen in all ionization states.

In Figure 14, we plot the oxygen deficit along with the lower limit mass estimate of Tumlinson et al. (2011) shown by the red dot-dashed line. At a stellar mass of 10^9.5 M_☉, our model showing the largest oxygen deficit (solid curve) is consistent with the lower limit estimates of Tumlinson et al. (2011). These independent results rule out all except the highest oxygen deficit model (solid line in Figure 14). A model with significant and ubiquitous outflows in star-forming galaxies is required for consistency with independent observations. As discussed in Section 8.4, our models suggest that in the case of significant outflows the oxygen deficit will scale (nearly) linearly with stellar mass. A linear scaling is consistent with the result found by Kirby et al. (2011) from their analysis of metal mass loss in Milky Way dwarf galaxies. At 10¹¹ M_☉, our model exceeds the lower limit estimate of Tumlinson et al. (2011) by ∼2 orders of magnitude.

Zoom In Zoom Out Reset image size

Gas propelled to high velocities may escape from galaxies all together. Estimates of the escape fraction of outflowing material have been difficult to determine accurately mainly due to lack of constraints on halo drag (Veilleux et al. 2005). It may be that detailed estimates of the oxygen mass in the CGM of galaxies will account for the ∼2 orders of magnitude greater oxygen deficit at M_* ∼ 10¹¹ M_☉ implied by our models as compared to lower limit estimates of Tumlinson et al. (2011). However, if the oxygen deficit is not resolved by oxygen found in the halos of galaxies, a natural reservoir for the remainder of the oxygen deficit would be the intergalactic medium (IGM). This scenario is consistent with theoretical models where large-scale galactic outflows transport metals from the galaxy to the IGM (e.g., Springel & Hernquist 2003; Oppenheimer & Davé 2006; Kobayashi et al. 2007; Oppenheimer et al. 2011). However, the kinematics of oxygen in the CGM of local star-forming galaxies suggest that most of the oxygen is gravitationally bound (Tumlinson et al. 2011).

A common method of estimating the escape fraction is to compare the outflow velocity to the escape velocity of the gravitational potential well of the galaxy. Early theoretical work suggested that mass loss occurs only in galaxies with shallow potential wells (e.g., Larson 1974; Dekel & Silk 1986; Mac Low & Ferrara 1999; Ferrara & Tolstoy 2000). However, recent observational and theoretical work has suggested that mass loss driven by AGN or star formation occurs in significantly higher mass galaxies as well (Strickland et al. 2004; Murray et al. 2005; Weiner et al. 2009; Steidel et al. 2010; Fischer et al. 2010; Feruglio et al. 2010; Sturm et al. 2011; Aalto et al. 2012).

It is not clear whether low-mass or high-mass galaxies are responsible for the IGM enrichment. Though low-mass galaxies may lose a larger fraction of their metals, they may not be a significant source of IGM enrichment due to their low rates of metal production (Martin et al. 2002; Kirby et al. 2011). On the other hand, due to the deep potential wells of large galaxies, metals may not be efficiently ejected into the IGM. If the latter is true, our model would predict that the CGM and/or outer disks of massive star-forming galaxies (10¹¹ M_☉) should contain substantially greater quantities of oxygen than their lower mass counterparts (10⁹ M_☉). Surveys of the CGM and studies of the outer disks of star-forming galaxies currently under way will provide important clues to resolving the oxygen deficit implied by our models.

One possible interpretation of the oxygen deficit is that it has been driven out of the galaxy ISM by galactic outflows. If oxygen is expelled via outflows, we expect an even larger mass of gas expelled since outflowing gas will not be pure oxygen. We can estimate the total baryonic mass loss in two ways. First, we adopt an enriched wind model assuming that the metallicity of outflowing gas is equal to the Type II SN oxygen yield (see Section 8.2). We adopt a value for P_o = 0.007, which is ∼1.6 Z_☉. This value is consistent with estimates of Martin et al. (2002) for the metallicity of a starburst-driven wind in NGC 1569. Assuming that outflows are enriched, we can estimate the total baryonic mass loss by dividing the oxygen deficit by the oxygen yield. This is given by

$$\begin{equation} \Delta M_b = \frac{M_T^o - M_g^o - M_\ast ^o}{P_o} = \frac{\Delta M_o}{P_o}. \end{equation} \tag{ 29 }$$

ΔM_b is the total amount of baryonic mass loss, and ΔM_o is the oxygen deficit.

We also derive an upper limit estimate of the total baryonic mass loss by assuming a uniform wind model. In the uniform wind model, the metallicity of the outflowing gas is the same as the ISM. In general, the metallicity of the outflowing gas is lower than the yield, and therefore more gas is required to outflow in this model in order to yield the oxygen deficit derived in Section 8.4. In the uniform wind model estimate P_o in Equation (29) is replaced by an SFR-weighted metallicity, Z_ψ. The SFR-weighted metallicity and the total baryon loss are given by

$$\begin{equation} Z_\psi = \frac{\int _{t_i}^{t} \Psi Z_g dt}{\int _{t_i}^{t} \Psi dt} \end{equation} \tag{ 30 }$$

and

$$\begin{equation} \Delta M_{b\psi } = \frac{\Delta M_o}{Z_\psi }, \end{equation} \tag{ 31 }$$

and ΔM_bψ > ΔM_b.

In the top panel of Figure 15, we plot the total baryonic mass loss as a function of stellar mass for galaxies in the local universe. The black and red curves are estimates adopting an enriched outflow (Equation (29)) and uniform wind model (Equation (31)), respectively. In the bottom panel of Figure 15, we plot the logarithm of δ_b, which we define as

$$\begin{equation} \delta _b = \frac{\Delta M_b}{M_g + M_\ast }. \end{equation} \tag{ 32 }$$

δ_b is the ratio of the total baryonic mass loss to the baryonic mass of the galaxy, where the baryonic mass of the galaxy is given as the sum of the gas and stellar mass. δ_b represents the ratio of the mass of gas to be cycled in and out of the galaxy compared to the current baryonic mass of the galaxy.

Zoom In Zoom Out Reset image size

In the currently accepted cosmological model (ΛCDM), the universal baryon fraction is precisely determined from the cosmic microwave background, the observed baryon acoustic oscillations, and the Hubble constant. The universal baryon and dark matter densities revealed by these observations are given by Ω_b = 0.0456 ± 0.0016 and Ω_c = 0.227 ± 0.014, respectively (Komatsu et al. 2011). Studies attempting to account for the baryons find that only a fraction of the expected baryons are observed in the low-redshift universe (Fukugita et al. 1998; Fukugita & Peebles 2004; Nicastro et al. 2005; Sommer-Larsen 2006; Shull et al. 2011).

Only about a tenth of the baryons are found in the stars and gas of galaxies (Bell et al. 2003a). While the Lyα forest at low redshifts can account for another ∼30% (Penton et al. 2004; Sembach et al. 2004), the majority of baryons are still missing. Some cosmological simulations favor the WHIM as the repository of the missing baryons (e.g., Cen & Ostriker 1999, 2006; Davé et al. 2001; Oppenheimer et al. 2011). However, observations of the hot gas at 10⁵–10⁷ K composing the WHIM remain tentative, and no compelling evidence for the detection of this phase yet exists (Bregman 2007). The distribution of the WHIM material is unknown, and hot halos of massive galaxies are considered as a possible reservoir for substantial fraction of the missing baryons (Cen & Ostriker 2006; Tang et al. 2009; Kim et al. 2009), though this remains controversial (Anderson & Bregman 2010).

An open question is what fraction of the missing baryons were ejected from galaxies through feedback processes and what fraction never accreted in the first place. We can compare our estimates of the total baryonic mass that can be associated with galaxies (baryonic mass plus the baryonic mass loss) with the expected baryon content of galaxies and their halos inferred from cosmological estimates. The stellar-to-halo mass (SHM) relation parameterizes the relationship between stellar mass of galaxies and the dark matter halos in which they reside. Moster et al. (2010) develop a statistical approach whereby halos and subhalos are populated in N-body simulations with the requirement that the observed stellar mass function be reproduced. They parameterize the SHM as

$$\begin{equation} \frac{M_\ast }{M_h} = 2 \left(\frac{M_\ast }{M_h}\right)_0 \left[ \frac{M_h}{M_1}^{-\beta } + \frac{M_h}{M_1}^{-\gamma } \right]^{-1}. \end{equation} \tag{ 33 }$$

Here M_h is the halo mass and (M_*/M_h)₀, M₁, β, and γ are free parameters. The relation evolves with redshift, and the parameters are given by

$$\begin{eqnarray*} \mathrm{log} M_1 (z) & =& 1.07 \cdot (1+z)^{0.019} \nonumber \\ \left(\frac{M_\ast }{M_h}\right)_0 (z) &=& 0.0282 \cdot (1+z)^{-0.72} \nonumber \\ \gamma (z) &=& 0.556 \cdot (1+z)^{-0.26} \nonumber \\ \end{eqnarray*}$$

and

$$\begin{equation} \beta (z) = 1.06 + 0.17 z. \end{equation} \tag{ 34 }$$

We can estimate the expected baryon content of galaxies from the universal ratio of baryonic to dark matter given by f_bc = Ω_b/Ω_c = 0.201 ± 0.014.

Figure 16 demonstrates the missing baryon problem for galaxies. The dashed cyan line is the current total baryonic (gas + stellar) mass content of galaxies plotted as a function of stellar mass. The dot-dashed blue line is the expected total baryonic mass content of galaxies assuming the universal baryonic-to-dark-matter ratio and the SHM relation of Moster et al. (2010) and is given by M_b = f_bc M_h. The solid black and red curves are the sum of the baryonic mass and the baryonic mass loss using the two estimates given by Equations (29) and (31), respectively. The solid black and red curves can be interpreted as estimates of the baryon content associated with galaxies and are the amount of baryons currently found in local star-forming galaxies plus what was once in the galaxies but has since been cycled out through outflows.

Zoom In Zoom Out Reset image size

The baryon content associated with galaxies is still substantially lower than the inferred content from cosmology. The straightforward interpretation of Figure 16 is that the missing baryons were never accreted onto local star-forming galaxies, and unless a large reservoir of gas is found in the halos of star-forming galaxies, it is likely that the missing baryons reside in the IGM. Using observationally motivated constraints for the mass and radii of hot halos, Anderson & Bregman (2010) come to a similar conclusion.

From our estimate of the baryon mass loss we can derive an effective mass loading factor, η, which is given by

$$\begin{equation} \eta = \frac{\Delta M_b}{M_\ast /(1-R)}. \end{equation} \tag{ 35 }$$

The quantity ΔM_b is the total amount of baryonic mass loss, and M_*/(1 − R) represents the total amount of star formation. In the literature the instantaneous mass loading factor is defined as $\eta _i = \dot{M}_w/\dot{M_\ast }$. The instantaneous mass loading factor is the mass loss rate divided by the SFR. The effective mass loading as we have defined it in Equation (35) is the instantaneous mass loading factor, η_i, averaged over the star formation history. In Figure 17, we plot the effective mass loading factor as a function of stellar mass. Here, we have adopted R = 0.5. Assuming a smaller value for R would lower the estimate. Even adopting our upper limit estimates of baryon mass loss, the effective mass loading for star-forming galaxies in the local universe is <1.

Zoom In Zoom Out Reset image size

A complete and self-consistent theory of galaxy evolution will require a detailed account of galaxy growth and chemical evolution along with physical mechanisms governing these processes. In this contribution, we present empirical models attempting to self-consistently integrate chemical evolution and galaxy growth. Our self-consistent approach is complementary to cosmological simulations and semi-analytical models, and the self-consistent census approach presented here should be compared with those models. Both theoretical and observational advances are crucial to constraining the model results developed in this study. The method of analysis used in this study provides useful tests for consistency of a diverse set of observations and theories related to chemical properties of galaxies. Here we address future prospects for improvement.

One of the greatest outstanding astrophysical problems is the large uncertainty in the absolute calibration of the nebular abundance scale. High- S/N observations of a large sample of H ii regions covering a broad range of physical parameters will be extremely important in statistically establishing and testing diagnostics. Such observations will be crucial in developing empirical calibrations relying on recombination lines, which are thought to be less susceptible to effects of temperature variations associated with auroral lines (Esteban et al. 2002; Peimbert & Peimbert 2005; Bresolin 2007). A complementary approach will be to use the recently developed wide-field integral field spectrographs to observe nearby H ii regions in order to understand the discrepancy between empirically and theoretically calibrated strong-line methods. A well-calibrated diagnostic applied to large data sets investigating the MZ relation at higher redshifts and to lower stellar masses will provide important constraints for our understanding of chemical evolution and our census of oxygen in the stellar and gas phase.

Theoretical models incorporating stellar rotation and mass loss into computations of nucleosynthesis for stars at all masses are yet to be developed (Romano et al. 2010). Current models of stellar nucleosynthesis are not yet able to explain the full diversity of chemical abundance patterns observed. While the oxygen yields are much better constrained by models owing to its primary origin, there is still a factor of ∼2 discrepancy. A single, self-consistent model of nucleosynthesis able to reproduce the abundance ratios observed in the Milky Way and other nearby galaxies will likely alleviate some of the discrepancy, thus constraining the total oxygen production in star-forming galaxies.

Part of the uncertainty in chemical evolution models rests on the assumption of a particular IMF (Romano et al. 2005). Adopting a constant IMF is problematic for galaxies, where even a constant stellar IMF could lead to variations in the integrated galactic IMF (Kroupa & Weidner 2003), and studies of the IMF from integrated measurements of star-forming galaxies indirectly indicate variations (Hoversten & Glazebrook 2008; Lee et al. 2009; Meurer et al. 2009). A varying IMF would have important implications for chemical evolution models (e.g., Romano et al. 2005; Köppen et al. 2007; Calura et al. 2010). Despite the mounting evidence, no direct evidence of IMF variations in star-forming galaxies is currently available, leading to a lack of consensus on the constancy of the IMF. Distinguishing whether a universal IMF or IGIMF formulation provides a better description of large-scale star formation will be important for chemical evolution studies of galaxies. Large statistical studies of the integrated properties of galaxies observed over cosmic time, in particular observations of the far-UV and far-IR properties, which are now becoming available, will be useful in establishing any systematic variations of the IMF. In our study, tighter constraints on the IMF will alleviate much of the uncertainty in the total amount of oxygen produced and the return fraction.

Total gas masses are required to measure the absolute content of metals within galaxies from observed relative abundances. In the local universe, large surveys of atomic and molecular gas have been conducted (Helfer et al. 2003; Walter et al. 2008; Leroy et al. 2009; Saintonge et al. 2011). Molecular gas at higher redshifts has also been detected in star-forming galaxies (Daddi et al. 2008, 2010; Tacconi et al. 2010). The next generation of radio and submillimeter observatories such as the Atacama Large Millimeter/submillimeter Array and the Square Kilometer Array will revolutionize the study of cold gas in the universe, allowing us to probe larger samples to far greater depths. Currently, we are only able to estimate the total oxygen in the gas phase of local star-forming galaxies. Measurements of cold gas in nearby and distant star-forming galaxies taken together with a well-calibrated gas-phase metallicity diagnostic will allow us to track the total mass, not just relative abundance, of oxygen in the gas phase of star-forming galaxies over cosmic time.

The structure, content, and chemical composition of outflowing gas are highly uncertain. Understanding the physical properties of galactic winds is crucial for estimating the total mass of gas outflowing from galaxies. Observations reveal that galaxy scale outflows have a complicated multiphase structure (see Veilleux et al. 2005), and observations suggest that structure and composition vary in each of the phases (e.g., Tripp et al. 2011). Multi-wavelength observations are required. The cold gas component can be observed in absorption lines of neutral and low ionization state metals, and the current generation of integral field spectrographs on 8–10 m class telescopes and wide-field integral field spectrographs on 4 m class telescopes will provide census of galactic scale outflows in local galaxies and presence of galactic winds in the distant universe. Hotter phases can be observed using space-based UV spectrometers such as the COS on board Hubble. A crucial but inaccessible phase is the so-called wind fluid that drives the stellar winds. A hard X-ray telescope with high spatial resolution and sensitivity is required to observe this phase. For the foreseeable future, astronomers will likely have to rely on detailed theoretical models and indirect observations to understand the wind fluid.

COS will also provide important insight into the circum-galactic and warm-hot ionized mediums, both thought to be important repositories of baryons. Studies of these regions will begin to reveal the baryonic content of the hot halos and IGM surrounding galaxies. Observations of these regions along with a census of baryons associated with galaxies will be crucial to potentially identifying the large fraction of missing baryons in the local universe and resolving this long-standing problem (Bregman 2007). A benchmark for these studies is to resolve the oxygen deficit in star-forming galaxies presented here.

In this study, we have applied the best theoretical and observational constraints available in undertaking a census of oxygen in star-forming galaxies. We are unable to draw any strong quantitative conclusions from the models developed owing to the large uncertainties associated with both our adopted model parameters and observational inputs. Nonetheless, this study represents one of the first attempts to self-consistently account for the stellar mass growth and chemical evolution of galaxies. As theoretical and observational advances allow for ever-greater constraints, we hope that the approach taken in this study will prove to be an important ingredient in testing and developing a fully self-consistent theory of galaxy evolution.