A CRITICAL LOOK AT THE MASS–METALLICITY–STAR FORMATION RATE RELATION IN THE LOCAL UNIVERSE. I. AN IMPROVED ANALYSIS FRAMEWORK AND CONFOUNDING SYSTEMATICS

In this section, we describe the cuts we use to define the samples drawn from SDSS spectroscopic survey, and illustrate their effect on the properties of the sample in Figure 1. Definition of the sample in this paper follows the procedures adopted by M10 to select the star-forming galaxies. As in M10, we start with the MPA/JHU reduction of SDSS DR7 spectroscopic sample. MPA/JHU catalog¹¹ contains spectroscopic line measurements derived using a custom pipeline that yields continuum subtracted fluxes (T04). Full MPA/JHU SDSS DR7 sample consists of 928,000 galaxies (all sample numbers are rounded to the nearest thousand), but this number drops to 212,000 after the application of the M10 redshift (0.07 < z < 0.30) and Hα signal-to-noise ratio cuts (S/N(Hα) >25). Relatively high signal-to-noise cut on Hα flux was applied in order to yield usable signal in other, weaker, emission lines that are necessary for the measurement of metallicities, while the relatively high low-redshift limit was chosen to ensure that spectroscopic fibers cover large fraction of each galaxy (2 kpc at z = 0.07 and 7 kpc at z = 0.30). Line flux errors are taken as listed in the MPA/JHU catalog.

Zoom In Zoom Out Reset image size

Next we follow M10 and apply minor cuts to remove anomalously low ($F(\rm {H}\alpha)/F(\rm {H}\beta)<2.5$) and very high Balmer decrements, i.e., dust attenuations (A_V > 2.5), bringing the sample to 203,000. The purpose of these cuts is to eliminate unphysical or extreme dust attenuation corrections. We presume that M10 follow Nagao et al. (2006), who determine the attenuation at V band from the Balmer decrement using Cardelli et al. (1989) extinction law:

$$\begin{equation*} A_V = 7.23\log \frac{F(\rm {H}\alpha)/F(\rm {H}\beta)}{2.86}, \end{equation*}$$

with the attenuation affecting the Hα line being related to A_V by

$$\begin{equation*} A_{\rm {H}\alpha } = 0.818 A_V. \end{equation*}$$

Note that A_V is evaluated in the V band as a matter of convention, but it represents extinction in H ii regions, which is typically several times larger than the extinction of the stellar continuum (Calzetti et al. 2000; Charlot & Fall 2000). Galaxies where emission line flux is dominated by non-stellar emission are removed using the Kauffmann et al. (2003) criterion on the BPT diagram (Baldwin et al. 1981) line ratios, coupled with the log N2 <−0.2 cut (N2 = F([N ii]6584)/F(Hα)), leaving 160,000 galaxies.

In Figure 1, we show the effect of the cuts as reflected in the specific SFR versus stellar mass plane. Specific SFRs are derived from the UV–optical spectral energy distribution (SED) fitting and masses are taken from MPA/JHU catalog.¹² Figure 1(A) shows the initial sample. We see the well-known bimodality in star-formation properties, with mostly quiescent galaxies lying below log (SFR/M_*) = −11.5 and actively star-forming galaxies lying on a relatively narrow sequence (the "star-forming sequence" or the "main sequence"). The upper dotted line (repeated in all panels) represents the star-forming sequence fit from Salim et al. (2007, hereafter S07), and pertains to galaxies selected as star-forming according to Brinchmann et al. (2004, hereafter B04) criteria, and the lower line shows the SF sequence as determined in Salim & Lee (2012) using the Gaussian decomposition of the full sample into star-forming and passive galaxies. Application of the redshift cut (Figure 1(B)), applied following M10, removes the low-mass tail, especially among galaxies with less intense SF. As a result, for log  M_* < 10, the sample tends to retain galaxies above the mean SF sequence relation. Application of the remaining M10 selection criteria, most importantly the Hα S/N cut and the BPT selection, removes quiescent galaxies but also eliminates a significant fraction of galaxies on the SF sequence (Figure 1(C)).

In consideration of the potential biases described above, we have also constructed an augmented sample, which retains all the cuts as the M10 sample except that it allows redshifts extending down to 0.005 as long as the fiber contains at least 10% (following Yates et al. 2012) of the total stellar mass (fiber masses are available from the MPA/JHU catalog). Even at lowest redshift, the majority of galaxies pass the fiber mass fraction criterion, so we do not expect a bias due to incompleteness. Note that SDSS mass covering fraction is typically 30% (95% of galaxies have the covering fraction between 17% and 50%). The augmented sample facilitates fuller characterization of the trends at log M_* < 10 (plot not shown). The augmented sample contains 259,000 galaxies (a 60% increase).

In this work, we will use four metallicity indicators: (1) M10 metallicities—an average of R23 and N2 estimates, (2) Bayesian metallicities of T04, (3) N2O2, and (4) R23 with an O32 term. We will primarily use metallicity estimates determined following a procedure described in M10. We refer to these metallicities as M10 metallicities. M10 metallicities are the average of 12+log(O/H) estimates from two strong-line methods: R23 (the ratio of four oxygen lines ([O ii]3727 doublet and [O iii]4958, 5007) to Hβ, Pagel et al. 1979) and N2 (the ratio of [N ii]6584 to Hα, van Zee et al. 1997). R23 and N2 metallicities were derived using relations from Maiolino et al. (2008), which calibrate various line ratios with respect to Kewley & Dopita (2002) theoretical metallicities of SDSS galaxies. All emission lines were corrected for dust attenuation using Hα-to-Hβ Balmer decrement and Cardelli et al. (1989) extinction curve. Following M10, we allow galaxies with log R23 <0.9 and log N2 <−0.35 (the range of validity of Maiolino et al. 2008 calibrations) and remove those where the two metallicity estimates based on R23 and N2 differ by more than 0.25 dex. We confirm that neither of these cuts removes galaxies of any particular SFR range. The final M10-like sample (i.e., sample with M10 redshift cuts) consists of 141,000 galaxies (matching the number quoted in M10), while the final augmented sample (sample allowing z < 0.07) consists of 222,000 galaxies. The median redshift of the final M10-like sample is 0.11 and that of the final augmented sample is 0.08.

For comparison with previous work and to discuss potential biases affecting M10 metallicities, we will consider additional metallicity estimates. One is the metallicity estimate taken from the MPA/JHU catalog, which was derived following the methodology described in T04. T04 metallicities are different from most other estimates because they are not based on some line ratio calibrated against another metallicity estimate, but come directly from fitting the stellar-continuum subtracted spectra to emission-line photoionization model spectra containing multiple emission lines (four Balmer lines and eight forbidden lines). T04 metallicities are not publicly available for all galaxies for which we calculate M10 metallicities because the MPA/JHU catalog provides them only for galaxies classified as star-forming by B04 criteria. Those criteria required a S/N threshold of 3 in all four BPT lines, but with flux errors scaled up to account for the scatter in fluxes of repeat observations, thus effectively becoming equivalent to S/N ratio cuts of between 5 and 7. As a result, of 141,000 (222,000) galaxies in the final M10-like (augmented) sample 93,000 (163,000) have T04 metallicities (∼60%).

Recently, Juneau et al. (2014) revisited the analysis of repeat observations and showed that flux error scalings are much smaller for flux ratios (needed for BPT classification and metallicity determinations) than for absolute fluxes and are very close to one. Therefore, in this paper, we do not scale flux errors.

We will also derive metallicities using the N2O2 method (the ratio of [N ii]6584 to [O ii]3727) for which theoretical calibration is taken from Equation (7) in Kewley & Dopita (2002) and using the R23+O32 method (R23 method with an O32 ([O iii]4958, 5007 to [O ii]3727 line ratio) term McGaugh 1991), based on theoretical calibration from Equation 18 in Kobulnicky & Kewley (2004). Unlike the M10 metallicities, these two methods are expected either to be less dependent on the ionization parameter (N2O2) or to explicitly correct for it (R23+O32).

In this study, we will primarily use four measurements of SFR. Two pertain to SFR contained within the SDSS spectroscopic fiber, and two are integrated (total) SFRs. Following M10, we calculate fiber SFR from the Hα luminosity, dust corrected with a Balmer decrement and converted to SFR using Kennicutt (1998) conversion. Another fiber SFR estimate comes from the MPA/JHU catalog, determined according to B04 method of fitting photoionization models to six emission lines (including Hα and Hβ) simultaneously. As in the case of T04 metallicities, the method does not apply any explicit conversions to obtain the target parameter (metallicity or SFR), but performs Bayesian fitting of fluxes or flux ratios. B04 have shown that photoionization models imply that the conversion between the dust-corrected Hα and SFR is strongly metallicity (and therefore, indirectly, mass) dependent. As a result, B04 SFRs can be as much as 0.5 dex higher for the most metal-rich galaxies. Fiber SFRs, originating from nebular emission of very massive stars, measure instantaneous SFR (timescale <10 Myr).

First of the two total SFRs comes from the SED fitting of GALEX (Martin et al. 2005) ultraviolet (UV) and SDSS optical broadband photometry. We perform SED fitting on all SDSS spectroscopic sample galaxies covered by GALEX medium-deep imaging (50% of SDSS area). Fitting includes galaxies covered but not detected by GALEX, invariably non-SF galaxies, thus keeping the sample optically selected. SED fitting is performed using Bayesian approach by comparing the observed fluxes with a library of model fluxes obtained from stellar population synthesis models (Bruzual & Charlot 2003) using some priors on star-formation history and attenuated according to the Charlot & Fall (2000) dust extinction model. The methodology is described in detail in S07. Principal differences with respect to S07 are that we now use model priors optimized for star-forming galaxies (as described in da Cunha et al. 2008) and we produce model libraries more finely sampled in redshift (0.01 versus 0.05 in S07). SED SFRs, being mainly constrained by UV, are determined as averaged over the past 100 Myr.

Second type of total SFRs come from 22 μm (W4 channel) observations from WISE (Wright et al. 2010). We use the AllWISE catalog profile-fit W4 fluxes¹³ and extrapolate them using Dale & Helou (2002) IR SEDs calibrated with Marcillac et al. (2006) relations to obtain the total IR luminosity, which is then converted to SFR using the Kennicutt (1998) relation. Details can be found in Salim et al. (2009), where the same procedures were used on Spitzer 24 μm fluxes. Detections at 22 μm are available for 54% of the final sample, but this incompleteness does not produce strong biases in SSFR–M_* parameter space (Figure 1(D)). Mid-IR emission is typically assumed to originate from dust-enshrouded young populations that also give rise to the UV emission and therefore are expected to trace current (∼100 Myr) SFR. However, Salim et al. (2009) have shown that mid-IR may have significant contributions from older populations even in actively star-forming galaxies, so that the timescale over which mid-IR measures SF may be effectively on the order of few gigayears.

We will also be considering specific SFRs. Fiber SSFRs are obtained by normalizing fiber SFRs by fiber stellar masses from the JHU/MPA catalog. SED SSFRs are derived directly from the SED fitting and mid-IR SSFRs are normalized by the total stellar mass from the JHU/MPA catalog. As in M10, we use stellar masses from the JHU/MPA catalog, which were derived from SDSS photometry (following but independent from S07). These masses are in a very good agreement (scatter in difference ∼0.08 dex, with <0.01 dex overall offset) with the masses obtained from our SED fitting.

For all MPA/JHU measurements (metallicity, SFRs, stellar masses) and for parameters from the SED fitting, we use medians of the probability distribution function as fiducial parameter values.

SFRs and stellar masses were either derived with or converted to Chabrier IMF.

In this section, we set the stage for the remainder of our analysis by identifying the optimal and the most physically motivated parameter driving the secondary dependencies in the MZR. We then use this knowledge to propose an intuitive framework for investigating the M_*–Z–SFR relation in Section 3.3.

A useful way to think of the M_*–Z–SFR relation is that it is represents an extension of a familiar MZR. Thus, a commonly used way to illustrate the SFR dependence is to show average (or median) mass–metallicity tracks of galaxies selected to lie in different bins of absolute SFR (e.g., M10, Yates et al. 2012; Andrews & Martini 2013; Nakajima & Ouchi 2014). In this section, we will show that the absolute SFR is not the optimal quantity with which to characterize MZR residuals. The choice to use absolute SFR as a secondary parameter in the mass–metallicity relationship is probably motivated by the fact that the M_*–Z–SFR relation has originally been formulated by M10 and LL10 as the relation between the mass, metallicity, and SFR. Instead, we argue that this should rather be the relative SSFR—the difference between galaxy's SSFR and the SSFR typical for galaxies of that mass (i.e., the offset from the star-forming sequence).

To show this, we will contrast MZRs (Figure 2) of samples selected to have the extreme (highest 2.5% (top panels)) and the lowest 2.5% (lower panels)) values of (1) absolute SFR (left panels), (2) absolute SSFR (middle panel), and (3) relative SSFR (right panels) to the MZR of the general population of galaxies (repeated green band in each panel is the 90 percentile range of the metallicities of the overall M10-like sample). Median trends of top (bottom) samples is shown as a purple (red) line, and of the general population as a white line. In Section 3.1, we will only be using M10 metallicities and SFRs derived in SDSS spectroscopic fibers, as derived in M10. Thus our Figure 2 is directly comparable to Figure 1 in M10.

Zoom In Zoom Out Reset image size

We begin this exercise by exploring MZR offsets due to absolute SFRs. We see that the galaxies with the highest SFRs (Figure 2(A)) show no offset with respect to the overall sample at the highest masses (log M_* > 11), while at log M_* = 10 the offsets have increased to 0.2 dex, and stay that large at lower masses. Since SFRs on average increase with mass (e.g., B04), the highest SFRs are preferentially found among the more massive galaxies. Likewise, the galaxies with the lowest SFRs will be found in the region of lower masses (panel (B)). Galaxies with the lowest SFRs are systematically offset above the median overall MZR, by some 0.05 dex, regardless of the mass.

Given that the SFR to first order simply scales with mass, it is justified, and physically more motivated (e.g., Ellison et al. 2008; Yates et al. 2012), to explore whether the SFR normalized by mass, i.e., the specific SFR—an indication of the current SF activity with respect to the past average—would produce stronger offsets across the mass range than just the SFR. We again show M_*–Z plots, but now for 2.5% of the sample with highest and lowest SSFRs in the sample (Figures 2(C) and (D)). Since on average, the specific SFR declines with increasing mass (e.g., S07, also Figure 1), we now have the reverse situation that the highest SSFRs are found among the lower-mass galaxies and the lowest SSFRs are among the more massive. Median offsets for the intense star-formers remain large at high masses, but they are somewhat smaller than in the case of top SFRs for the lowest masses. For galaxies with low SSFRs (panel (D)), there is not much overlap in the mass regime with low SFRs (panel (B)), and the median MZR also sits some 0.05 dex above the overall median.

One can see that considering galaxies selected by either SFRs or SSFRs implicitly introduces a mass selection. A true secondary parameter driving offsets in MZR should apply to galaxies of all masses. Thus we now explore samples selected by the level of SSFR relative to what is typical at a given mass. Relative SSFR (Δlog SSFR) is defined as

$$\begin{equation*} \Delta \log {\rm SSFR} = \log {\rm SSFR} - \langle \log {\rm SSFR}\rangle _{M_*}, \end{equation*}$$

where $\langle \log {\rm SSFR}\rangle _{M_*}$ is the median or mean of log SSFR of galaxies having a mass M_*. In this work we use medians in 0.15 dex wide mass bins. Relative SSFR can be visualized as the offset of a galaxy from the star-forming sequence (vertical distance in Figure 1), a measure also called "SFR excess" (Schiminovich et al. 2007) or "starburstiness" (Elbaz et al. 2011). The relative SSFR is a "natural" parameter to characterize star-formation activity and is becoming increasingly used in the recent literature (e.g., Magdis et al. 2012; Woo et al. 2013). In Figures 2(E) and (F), we now show galaxies with the highest and lowest 2.5% relative SSFRs in each mass bin. These galaxies typically have SSFRS that are five times higher/lower than typical values at that mass |Δlog SSFR| > 0.7 dex). The quantity of points at different masses now reflects the mass distribution in the overall sample. We notice that some offset for highly star-forming galaxies (panel (E)) is now present even at the highest masses, which was not the case when absolute SFRs were considered (panel (A)). At lower masses, the offsets are at least as strong as they were in panel (A), but now the selection includes more galaxies. Similarly, for galaxies with the lowest relative SSFRs (panel (F)), the offsets are as large as in the case of either the lowest SFRs or SSFRs, but spanning the full range of masses.

From this section, we conclude that MZR offsets are more naturally characterized by the difference between the logarithm of galaxy's SSFR with respect to a typical log SSFR at a given mass, rather than the absolute SFRs. In hindsight, this may seem fairly obvious, but this point has not has not previously been clearly made. The interpretation of the MZR offsets in terms of the variations of relative SSFR is more in accordance with FMR's evolutionary sense: at higher redshifts the SF sequence appears to shift upward without much change in the slope (e.g., Speagle et al. 2014), so considering local samples at some distance (offset) from the SF sequence is more analogous to a high-redshift selection.

Does the above mean that the functional form of Z(M_*, SFR) should feature SSFR rather than SFR? In the case of a linear relationship and total measurements, the two formulations are formally equivalent. However, for fiber measurements, the difference is crucial. Fiber SFR, used in M10's formulation of FMR (their Equations (2), (4), and (5)), is not a meaningful physical quantity since it depends, to the zeroth order, on galaxy distance: on average, SDSS fibers cover 25% of SF at z = 0.07, but 65% at z = 0.30. On the other hand, fiber specific SFR (fiber SFR normalized by mass in the fiber) is a perfectly valid physical quantity, representing the intensity of SF in the same physical region of the galaxy in which the metallicity is measured and is therefore probably even preferred to the total specific SFR for the purposes of M_*–Z–SFR analysis. Future studies should report their fiber or slit SSFRs to allow for a more direct comparison with SDSS. Surprisingly, M10 find that their least-scatter projection, formulated with fiber SFRs, agrees with high-redshift measurements, even though fiber SFRs are distance-dependent and are on average 0.6 dex smaller than the total SFRs (but with a substantial, distance-dependent scatter). Altogether, it makes more sense to formally describe M_*–Z–SFR relations using specific SFRs. Indeed, the recent analytical model of Lilly et al. (2013, e.g., their Equation (40)) finds the metallicity to be the function the SSFR (see discussion in Section 5).

Guided by the inferences made in this section, in the rest of the paper we will consider the relative SSFR as the primary independent variable for galaxies of a fixed mass.

Previous work onM_*–Z–SFR relation in SDSS used exclusively SFRs derived based on the emission lines measured in spectroscopic fibers. As some of those same line measurements are involved in the determination of metallicity, this opens a concern that the M_*–Z–SFR relation may to some extent be the result of spurious correlations between the measurements of SFR and the metallicity. Therefore, in this paper, we will examine the relationship with two completely independent total SFRs, based on integrated fluxes. Furthermore, different SFR indicators are sensitive to SF over different timescales, which in principle may be more or less strongly tied to the changes in the metallicity.

In this section, we investigate four SFR indicators: two measured in fibers (both based on emission lines, but one following M10's common methodology of Balmer decrement-corrected Hα luminosity, and the other using B04's more sophisticated methodology of modeling simultaneously multiple emission lines) and two total measurements (one based on the UV/optical SED fitting, and the other based on 22 μm mid-IR luminosity from WISE). The results are presented in Figure 2 (right panels—M10 SFRs), and are continued in Figure 3 (left panels—B04 SFRs; middle panels—SED SFRs; right panels—mid-IR SFRs). In Figure 3, we continue to plot the galaxies with 2.5% highest/lowest relative SSFRs (top/bottom panels), and contrast them to the general population of galaxies (green band). We continue to use M10 metallicities.

Zoom In Zoom Out Reset image size

We begin with the comparison of M_*–Z plots based on M10 fiber SFRs (Figures 2(E) and (F)) to B04 fiber SFR (Figures 3(A) and (B)). Below log M_* = 10.3 there is not much difference in MZRs of either the high or the low SSFR samples selected by either M10 or B04 SSFRs. Above log M_* = 10.5, galaxies selected by B04 SSFR show no offset with respect to the overall metallicity. Since B04 SFRs are more sophisticated than M10 SFRs and therefore presumably more accurate, this suggests that there is very little or no SFR dependence in the MZR above log M_* = 10.5.

Next we look at SED SFRs—total SFRs determined from the UV/optical broadband fluxes, and primarily constrained by the UV. The MZR of the most intensely star-forming galaxies (Figure 3(C)) shows approximately two times smaller offset below log M_* = 10.3 with respect to equivalent relations based on either fiber SFRs. Above, log M_* = 10.5, the offset is nearly gone, as in the case of B04 fiber SFRs. For galaxies selected to have the lowest relative SSFRs (Figure 3(D)), the offset is likewise smaller than for fiber SFRs and basically not present above log M_* = 10.3, again similar to B04 SFRs. Before drawing any conclusions, we examine high/low samples selected by mid-IR SFRs. For intense star-formers (Figure 3(E)), the offset in MZR is, surprisingly, again as strong as it was for M10 fiber SFRs and is clearly present even at log M_* = 10.5, which was not the case for either B04 fiber SFRs or for the total SED SFRs. For galaxies furthest below the SF sequence (but detected at 22 μm), the offset is quite small (0.02 dex), as in the case of SED SFRs, and unlike M10 SFRs. We have verified that this small offset is not because of 22 μm detection limit potentially eliminating lowest star-formers—if we produce M10 SFR selected M_*–Z plot (not shown) but only for galaxies also detected at 22 μm, we still obtain an offset as large as the one in Figure 2(F).

To conclude, we find that the dependence of MZR on SSFR is definitely present using all SFR indicators, whether they be fiber or total, so the M_*–Z–SFR relation cannot be merely an artifact of correlated measurements. The degree of offsets vary, but at least one SSFR of both types (fiber and total) shows equally strong offsets for highly SF galaxies. It is interesting that the dependence is not significantly weakened when using total SSFRs, considering that the metallicity is measured in the fiber. Note that if an SSFR measure has a high uncertainty, it will not be able to accurately identify galaxies that are truly the highest star-formers and therefore the most discrepant in metallicity, which may explain why the trends are weaker using SED SFRs. This emphasizes the need to perform comparison of samples at different redshifts using indicators of similar accuracy, and, preferably, of the same type. Results further suggest that the fact that different indicators trace SF over different timescales does not appear to affect the trends significantly. Finally, we confirm previous findings (Ellison et al. 2008, M10) that the MZR offsets are more pronounced at lower masses, and that they are close to zero when approaching log M_* = 11.

In previous sections, we used familiar M_*–Z plots and contrasted samples of galaxies selected to be extremes in terms of star-forming properties. We now wish to include all of galaxies in the analysis. We have also seen that the MZR offsets are strongly mass dependent. We therefore need to analyze galaxies of different mass ranges separately. Following the framework set up in Section 3.1, we do so by plotting the metallicities against the relative SSFR. These Z–ΔSSFR plots for galaxies belonging to some mass bin embody our non-parametric methodology for exploring the M_*–Z–SFR relation.

In Figure 4, we show four Z–ΔSSFR plots, for one of 0.5 dex wide mass bins centered on log M_* = 9.5, 10.0, 10.5, and 11.0. Results do not depend strongly on the width of mass bins. We continue to use M10 metallicities and revert to M10 fiber SSFRs, but from now on use the augmented sample (one extending to lower redshift, z = 0.005, than what was used in M10), which allows for better characterization of the low-mass regime (Section 2.1). Relative SSFR for a given galaxy is determined as the difference of its SSFR with respect to the median SSFR in 0.15 dex wide mass bin. Positive relative SSFRs correspond to galaxies sitting above the star-forming sequence. Almost identical results would be obtained if, instead of using running medians, the relative SSFR was calculated as the difference between SSFR and the SSFR corresponding to the linear fit of log SSFR versus log  M_* at that mass (such as those in Figure 1). Figure 4 now allows us to explore dependence of Z on SSFR for all galaxies, not just those at the extreme of SSFR distribution.

Zoom In Zoom Out Reset image size

We start from the mass bin centered on log  M_* = 9.5 (Figure 4(A)), where the metallicity trend is stronger than in higher mass bins. Fitting the linear relationship (purple line) yields a slope κ = d(12 + log (O/H))/dlog (SFR/M_*) of −0.18. However, running medians reveal that the SFR dependence is considerably stronger above the SF sequence (a linear fit would yield, κ_high = −0.26), then in its core and below it (κ_low = −0.12). Moving on to the next mass bin (Figure 4(B)), we find that the overall dependence on SFR gets weaker (κ = −0.12), but again the slope is steeper above the SF sequence than in its core and below it. Remarkably, the slope above the SF sequence is as steep as in the lower mass bin (values of slopes are given in Table 1). At log M_* = 10.5 (Figure 4(C)), the overall slope is only κ = −0.05, and while it is again steeper above the SF sequence (and still as steep as at lower masses), the steepening does not begin until 0.7 dex above the sequence and consequently encompasses only a small number of highly star-forming galaxies. Finally, in log M_* = 11.0 bin (Figure 4(D)), we stop seeing different behavior above and within/below the SF sequence, with a rather shallow overall trend of κ = −0.03.

Are the results presented here valid for total SSFRs? We list the slopes of metallicities against the offset from the SF sequence based on mid-IR SFRs in Table 1). The results are remarkably similar to those based on fiber SFRs, with the most notable difference being that the slope in the highest mass bin is even weaker. This basically confirms the conclusion from Section 3.2 that the differences in timescales of SFR indicators and whether they pertain to fiber or integrated measurements do not lead to great differences in metallicity's dependence on SSFR.

To conclude, in this section, we have demonstrated that the metallicity's dependence on SSFR is stronger above the SF sequence (as remarked by M10). And while the dependence for the bulk of galaxies (those in and below the SF sequence) weakens as the mass goes up, the dependence for intense star-formers (lying at ≳ 0.6 dex from the SF sequence) stays very similar, suggesting that the different mechanisms drive the M_*–Z–SFR relation depending on SF intensity, which itself may be related to the existence of different modes of SF (e.g., quiescent versus merger-driven). Some studies, especially at z ≳ 1, distinguish populations with high relative SSFR (e.g., Δ log SSFR >0.3, Elbaz et al. 2011) as having a special mode of star formation associated with mergers and label them "starbursts." However, at z ∼ 2, the starbursts produce a clear excess above the Gaussian distribution of log SSFR and dominate already at Δ log SSFR =0.6 (Rodighiero et al. 2011). However, while we find the high-end distribution of log SSFR in SDSS to eventually depart from a Gaussian, this excess does not become dominant (twice as high as the Gaussian) until Δ log SSFR >1.0 dex (for M10 SFRs, at log M_* = 10), well above the onset of break in the Z–SSFR relation. Therefore, we refrain from equating the two-mode behavior in metallicity trends with the normal SF versus merger-driven starburst distinction at this time, but do not rule out such a connection.

Important implication of this finding is that describing or extrapolating the SFR dependence using simple linear trends (equivalent to assuming a flat "fundamental plane" in LL10 or a single preferred projection of the FMR in M10) could lead to inconsistencies when comparing to high-redshift samples, as the local trends will be dominated by galaxies that show weaker SFR dependence. We see that using a single linear trend can produce discrepancies as large as 0.2 dex in the metallicity for the most active star-formers. Such difference can give very different character to the interpretation of high-redshift data. Recently, Maier et al. (2014) showed (their Figure 5) that using different descriptions of the M_*–Z–SFR relation by M10 (second-order polynomial versus the plane of the least scatter) produce different extrapolations for high-SSFR samples. At a given mass, a better way to parameterize the SSFR dependence is with a broken linear fit—one for galaxies above the SF sequence and other for the rest. Ultimately, to test whether high-redshift samples follow the local M_*–Z–SFR relation (and to determine whether they exhibit the dependence on SSFR internally), it is best to show both the local and the high-redshift data on metallicity versus SSFR plots, within some mass bin (lower masses will provide stronger diagnostic as long as samples are not too small). We apply this methodology in S. Salim et al. (in preparation) to test whether local M_*–Z–SFR relation is consistent with z ∼ 2.3 measurements.

The existence of the (S)SFR-dependence of MZR implies that if the (S)SFR was accounted for, the scatter in the relation would decrease. This decrease was presented in M10 as being dramatic. For example, their Figure 5 shows that the scatter in metallicity becomes a factor of almost three smaller (goes from 0.055 dex to 0.02 dex). However, what is perhaps not sufficiently appreciated is that this calculation in M10 pertained to the reduction of the rms residuals of median-binned values of metallicity around the best-fitting surface, and not of the individual galaxies. Different studies have claimed conflicting results regarding the scatter of individual galaxies (non-binned samples). Ellison et al. (2008) reported that SSFR is not an important cause of scatter in MZR (a 10% reduction) and more recently, Pérez-Montero et al. (2013) quoted a 0.01 dex reduction. These results are at odds with the M10 stated reduction from 0.08 dex to 0.05–0.06 dex.

The methodology applied in Section 3.3 allows us to directly address the question of the reduction of scatter. Figure 4 shows that the scatter in metallicity is high even at a fixed mass and fixed SSFR (metallicity axis in Figure 4 spans the exact same range as in M_* − Z plots: Figures 2 and 3). Even at low masses (log M_* = 9.5) where the SFR dependence is the strongest, the overall scatter (standard deviation) in metallicities of 0.12 dex is reduced only to 0.10 dex (scatter around the median), a 20% reduction. The effect is even more modest at higher masses, where the dependence on SFR is weaker. The values of the standard deviation of metallicities before and after accounting for the SFR are given in Table 1. The uncertainty in the measurement of SSFR (∼0.2 dex) adds to the scatter in Z, but this contribution (0.2κ) does not account for more than 10% of the residual scatter. Based on this, we conclude that presenting the MZR using a quantity that combines the mass and SFR (as in M10 Figure 5) should not be expected to reveal conspicuous reduction in the scatter for typical samples of galaxies. On the contrary, the reduction is indeed very modest. This also means that the M_*–Z–SFR relation, as a three-dimensional "surface," is not particularly thin, despite common wisdom (e.g., Steidel et al. 2014), possibly inspired by M10's use of the binned points (e.g., their Figure 2).¹⁴

Finally, we turn to another question regarding the scatter in Z. T04 and Zahid et al. (2012) have previously noted that the metallicity scatter in MZR, which becomes quite large at lower masses, is not due to, for example, larger measurement errors, but is intrinsic. The question is whether it can be fully explained by dependence on SFR. As we have seen, taking into account the SSFR dependence reduces the scatter at log M_* = 9.5 from 0.12 only to 0.10 dex, which is still significantly larger than the typical metallicity errors at that mass (0.04 dex, T04). This suggests that some other parameter(s) may ultimately be more closely connected with metallicity regulation than the SFR.

In this section, we explore how M_*–Z–SFR relation is affected by the use of alternative metallicity indicators, primarily those based on T04 method. Since T04 metallicities that are available in MPA/JHU catalog also involve additional selection criteria, we first need to explore if such criteria alone affect the dependence on SFR. In the analysis, we will continue to use Z–Δ SSFR plots split by mass bins.

Given all the practical and theoretical difficulties concerning the measurement of metallicities (e.g., Andrews & Martini 2013 and references therein), it is of critical importance to understand if different methods of deriving them affect the conclusions regarding the existence and the character of M_*–Z–SFR relation. In all of our analyses up to this point, we have used metallicities derived according to the method of M10, i.e., an average of estimates based on R23 and N2 calibrations of Maiolino et al. (2008). T04 metallicities, being available as part of the MPA/JHU catalog, are very often used in the analysis of SDSS samples and have been used for the characterization of M_*–Z–SFR relation (Lara-López et al. 2010; Yates et al. 2012). However, T04 metallicities are not available for the entire sample that we considered so far (which was selected by requiring S/N(Hα) > 25). Instead, T04 metallicities are available for galaxies satisfying the condition of having relatively high S/N ratios in all four BPT lines (S/N ratio >7.3, 5.5, 4.5, and 6.0 in Hα, Hβ, [O iii]5007, and [N ii]6584, respectively), as explained in Section 2.2, which is fulfilled for two-thirds of our sample. To see whether multiple S/N ratio cuts alone bias the trends, we repeat Figure 4, i.e., we still use M10 metallicities, but now restricted to those galaxies for which T04 metallicities exist. The resulting trends, in same mass bins as before, are shown in Figure 5. Comparison with Figure 4 reveals that while the overall trends are similar, there are important differences. Namely, the median metallicities below the SF sequence (Δ log SSFR <0) are now lower by up to 0.05 dex. There is no change above the SF sequence. In other words, the selection based on four emission lines leads to the preferential removal of galaxies with lower SSFRs and higher metallicities. The overall result is that the bulk trends (purple lines in Figure 5) become slightly weaker.

Zoom In Zoom Out Reset image size

Keeping in mind the biases introduced by S/N ratio cuts present in T04 sample, we now look at the SSFR trends using the actual T04 metallicities (Figure 6). The results are remarkably different compared to equivalent plots made with M10 metallicities (Figure 5). In the lowest mass bin (panel (A)), where the dependence on SFR was the strongest, it is now much weaker, with the overall slope of only κ = −0.05. More importantly, the sense of the Z(M_*, SFR) (that the more active galaxies at a given mass should have smaller metallicities) is only observed above the star-forming sequence (Δ log SSFR >0.5, and therefore involves a smaller fraction of galaxies. Bulk of the galaxies (those within the core of the star-forming sequence and below it) show no trend at all. Similar situation persists in higher mass bins. Anti-correlation between metallicity and SFR is observed only at ≳ 0.4 dex above the SF sequence, but for the bulk of the galaxies it is not present, and even turns into a weak positive correlation. The end result of using T04 metallicities is that the trends are no longer monotonic, which is inconsistent with the possibility that the M_*–Z–SFR relation is redshift invariant.¹⁵

Zoom In Zoom Out Reset image size

To disentangle the effects of T04 sample selection that requires high S/N ratios in multiple lines from those that are intrinsic to the T04 metallicity method, we have also calculated the metallicities, using the exact methodology of T04, for galaxies that do not pass T04 cuts. Thus we arrive at the sample that is selected with only the S/N ratio cut on Hα line. Investigating the metallicity trends with such sample we still find (plots not shown) that T04 metallicities lead to much weaker trends against SSFR than M10 metallicities, but not so much to produce a drop at low SSFRs that would lead to non-monotonic behavior. This, apparently, is due to the additional bias of having S/N ratio selections in multiple lines.

Differences between M10 and T04 metallicities in the context of M_*–Z–SFR relation were previously noted by Yates et al. (2012), who ultimately preferred the T04 metallicities. Considering that most studies determine metallicities using simple methods similar to those used in M10 and not using the complicated Bayesian fitting of full emission line spectrum as in T04, it is important to understand if the simple method of M10 is at fault.

One concern regarding the M10 method is that it is based on semi-empirical calibrations (Maiolino et al. 2008). Maiolino et al. (2008) calibrate various individual line ratios of a sample of SDSS galaxies against the metallicities obtained from the theoretical calibrations of Kewley & Dopita (2002). Thus, such calibration will follow the theoretical models on average, but not necessarily for parts of the sample that have properties different from the average. Specifically, their calibrations may not be valid for galaxies that have ionization parameters that are very different from what is typical at a given line ratio value. To test the possibility that the M_*–Z–SFR relation based on M10 metallicities is affected by biases due to the variations of ionization parameter, we also calculate metallicities using the N2O2 method, which is the least sensitive to ionization parameter of all simple methods (methods that employ one ratio of lines; Kewley & Dopita 2002). We show the N2O2 metallicity versus SF trends in Figure 7(A), but now only for one mass bin (log M_* = 10.5). The general sense of the trend is the same as it was with M10 metallicities: there is an anti-correlation for galaxies within the SF sequence, which becomes stronger above it. The overall slope is even slightly steeper than it was with M10 metallicities. There is no sign of flat or positively correlated trends as with T04 metallicities.

Zoom In Zoom Out Reset image size

In their choice to adopt T04 metallicities as more fiducial, Yates et al. (2012) put forward an argument that any method for deriving metallicities involving nitrogen (such as M10's, which is an average of N2 and R23) is possibly affected by saturation. To address this concern, we now consider an oxygen-based (R23) calibration of Kobulnicky & Kewley (2004), which at the same time accounts for differences in the ionization parameter (through its dependence on O32). The results (again for log M_* = 10.5 bin) are shown in Figure 7(B). While the median trend within the SF sequence is now somewhat weaker than for either the M10 or N2O2 metallicities, they do not go away as in the case of T04 metallicities. In other mass bins (plots not shown), both the N2O2 and the R23+O32 trends are again closer to trends using M10 metallicities than those using T04 metallicities.

Based on the analysis presented in this section, we conclude that the use of T04 metallicities, even when accounting for the biases due to having multiple-line selections, produces SFR dependencies that are much smaller than those produced using conventional line ratio methods. It should be kept in mind that even among the methods that yield similar results, the strength of the dependence on (S)SFR will differ depending on the metallicity method and calibration (Andrews & Martini 2013). This reiterates the point that the comparisons of different samples must be based on same metallicity method and calibration (in addition to using comparable (S)SFR estimates, and allowing for non-linear trends, as emphasized in previous sections).

Furthermore, in this section we find that the dependence on SFR will depend on the S/N ratio cuts applied to the lines as discussed in de los Reyes et al. (2014). This is especially true for galaxies with lower (S)SFRs. This can be avoided by using a S/N ratio cut only on one hydrogen line (which can be made sufficiently high so that other lines are also well measured).

The work that originally drew attention to the dependencies of MZR on other parameters, Ellison et al. (2008), claimed to have found not only MZR offsets due to SSFR, but also even stronger offsets resulting from galaxy physical sizes, in the sense that at a given mass larger galaxies have on average smaller metallicities. However, subsequent studies focused exclusively on the SFR aspect, neglecting the question of galaxy sizes, perhaps because the latter phenomenon had less clear intuitive interpretation, and because any result involving galaxy sizes in a sample where the physical scale and the fiber covering fraction span such wide ranges appears suspect. Furthermore, as pointed out by Ellison et al. (2008), simultaneous presence of both dependencies leads to the predictions that high-redshift MZRs should have much smaller evolution than what is observed, so the existence of a strong size dependence would potentially conflict with the basic idea of FMR—that it can account for MZR evolution.

Here we revisit the question of MZR's dependence on galaxy size by applying the methodology introduced in Section 3 (metallicity trends in 0.5 dex wide mass bins), but substituting SSFR with galaxy size. We use the same data for galaxy sizes (Simard et al. 2011) as used by Ellison et al. (2008), except that we take semimajor axis half-light radii in the g band (as opposed to r), based on single Sersic fits. The results for our augmented sample, and using M10 metallicities, are shown in Figure 8. Significant trends are seen in all mass bins, but they are stronger at lower masses. Bulk trends (purple fits) are not as strong as those versus SSFR, although they become more comparable at higher masses. The main difference is that the dependence versus SSFR has a two-mode behavior such that the trends above the SF sequence are much stronger, which results in the full range of median metallicities that is many times greater than the range of median metallicities due to the size dependence. For example, at log M_* = 9.5 the total span of metallicities due to the size is ∼0.1 dex, while it is ∼0.5 dex with respect to the SSFR (Figure 4(A)). Interestingly, the trend versus size reverses for very large galaxies (half-light size >10 kpc), but there are very few such galaxies compared to those of smaller size.

Zoom In Zoom Out Reset image size

One might have a concern that the size dependence is a just consequence of there being the primary dependence with SSFR. This would be the case if the size and SSFR were themselves strongly correlated in the sense that larger galaxies have higher SSFRs. However, this is not the case. We find (plots not shown) that the galaxies with higher SSFR on average have smaller sizes than galaxies with lower SSFRs (up to a factor of two). Furthermore, the dependence of metallicity on SSFR persists with the same intensity even when the galaxies are selected to lie in a very small range of physical sizes (plots not shown). Therefore, the two phenomena are independent. Finally, we confirm that the trends of metallicity with size remain when the sample is restricted to narrow redshift ranges (0.01), which demonstrates that the observed dependence is not an artifact of redshift-dependent covering fraction of the fibers.

The conclusion of this section is that the metallicity shows dependence on galaxy size that is independent from the dependence on SSFR and is not an artifact of aperture effects, nor of metallicity gradients within the galaxy (for the latter, see Ellison et al. 2008). However, the total trends are weaker compared to those with respect to SSFR, so one does not expect that the size dependence alone will produce much evolution in the MZR and thus affect the potential invariance of M_*–Z–SFR relation. Nevertheless, the size does appear to be a genuine contributing source of the dispersion (albeit small) in the MZR and is therefore in need of being tested with theoretical models.

The character of the M_*–Z–SFR relation, that at a given mass there is an anti-correlation between SFR and metallicity, has been brought into question by Yates et al. (2012; and to some extent Lara-López et al. 2013), who find that this dependence reverses above log M_* ≈ 10.2, such that the higher (S)SFRs are associated with, on average, higher metallicities. Such result, if correct, would basically preclude the "fundamental" aspect of the FMR. Namely, as we go to higher redshifts, and SFRs at a given mass rise, one would expect, based on Yates et al. (2012) "reversal," that the metallicities of low mass galaxies would be offset lower compared to local galaxies (as observed) and that the metallicities of higher mass galaxies should on average be located above the local MZR. No such evolution of MZR has been reported in the current literature. So either the local relationship between Z, M_*, and SFR does not at all hold at other redshifts (is not fundamental), or there is a problem with the finding that there exists a reversal of trends.

Yates et al. (2012) found the reversal using T04 metallicities. They have also considered, but eventually decided not to trust, the metallicities determined according to the M10 method, which, even in their analysis, did not show any evidence of the reversal (their Figure 1). Our analysis (Section 3.5) has shown that T04 metallicities intrinsically shows much weaker trends than other metallicity indicators, which are additionally exacerbated by T04 metallicities being available only for galaxies that fulfill S/N ratio criteria on multiple lines, leading to trends that are strongly non-monotonic, especially at higher masses (Figures 6(C) and (D)). Can this explain the results of Yates et al. (2012)? In Figure 9(A), we show metallicity versus SSFR for mass bin centered at log M_* = 10.5, where, according to Yates et al. (2012), the reversal should already take place. The sample selection used in this figure follows that of Yates et al. (2012), applying S/N > 5 cuts in Hα, Hβ, and [N ii]6584, redshift range of 0.005 to 0.25, and the requirement that fiber captures at least 10% of the r-band flux. Yates et al. (2012) use aperture-corrected total SFR from B04.¹⁶ Indeed, the general sense of the trend in Figure 9(A) (purple line) is one of a correlation (and therefore the reversal compared to the anti-correlation at lower masses). This is similar to what we have already seen when we discussed T04 metallicities (Figure 6), but with even fewer galaxies participating in the anti-correlation part of the trend, perhaps due to a different type of SFR or somewhat different sample cuts.

Zoom In Zoom Out Reset image size

One could argue for reversal being real by noting that a similar effect may exist in trends of dust extinction versus SFR. Namely, Zahid et al. (2013) find that dust extinction A_V (estimated from the Balmer decrement) in SDSS shows an anti-correlation with SFR at lower masses, which turns into a positive correlation above log M_* > 10.2. These trends are relatively weak and have a substantial scatter (∼0.5 mag). For SFR, Zahid et al. (2013) use B04 total SFRs. Interestingly, we confirm Zahid et al. (2013) results, but only when using total SSFRs. When we instead look at A_V versus fiber SSFR (either derived as in M10 or from B04), the mean trends have the same character irrespective of the mass: dust extinction rises with SSFR, reaches a peak, and then turns down above the SF sequence. The rising part of the trend becomes steeper with mass. Given that the Balmer decrement is determined in the fiber, it is more appropriate to compare it to the SSFR also measured in the fiber. Thus the results of Zahid et al. (2013) may not hold when the dust and the SSFR are measured in matching regions, possibly because the dust extinction and the SSFR do not scale alike.

T04 metallicities were also used in LL10, the other of the two papers that first reported the relationship between mass, SFR and metallicity. Like Yates et al. (2012), LL10 used total SFRs from B04, but also a much more restricted redshift range (0.04 < z < 0.10) and very high (>8) S/N ratio cuts in eight emission lines (four BPT lines, plus [O ii]3726, 3729 and [S ii]6717, 6731). We replicate LL10 selection and show the resulting metallicity versus SSFR plot at log M_* = 10.5 in Figure 9(B). The apparent correlation ("reversal") is now even stronger than in Yates et al. (2012). This is the result of lower redshift range which preferentially eliminates less frequent high-SFR galaxies that would have added weight to the anti-correlation, and due to the high S/N ratio cuts that tend to preferentially remove low-SFR galaxies with high metallicities. LL10 do not discuss the reversal in their original paper, but do confirm it subsequently (Lara-López et al. 2013). However, it must be noted that the high-mass reversal is entirely at odds with the concept of the "fundamental plane" introduced by LL10. The sense of the fundamental plane is always one of anti-correlation. This can be seen in Figure 9(B), where we plot the locus of LL10 fundamental plane as the orange line (it is also evident in Figure 13 (top left panel) in Lara-López et al. 2013). Furthermore, the fundamental plane requires the slope κ to be mass-independent (e.g., Figure 13 (top left panel) in Lara-López et al. 2013), which is obviously not the case no matter which metallicity or SFR indicator is used. The fact that LL10 were able to derive the fundamental plane despite using T04 metallicities that cause the apparent reversal, is because the plane was constrained by more numerous lower mass galaxies that dominate in their sample and for which the general trend is that of anti-correlation, even using T04 metallicities (Figure 6(A)).

The conclusion is that the apparent reversal in metallicity versus (S)SFR, as seen in Yates et al. (2012), is primarily the result of using the available T04 metallicities, and would not be seen using other metallicity estimates (or even with T04 metallicities if they were available for a sample not biased by multiple-line S/N ratio cuts). Furthermore, the reversal would have conflicted with the fundamental aspect of the local Z–M_*–SFR relation because it would no longer be able to explain the MZR evolution.

One potentially serious limitation of all studies of the relationships between metallicity, mass, and SFR that are based on SDSS data is that the metallicity measurements come from fiber spectroscopy, which covers part of the galaxy in a way that is redshift and galaxy-size dependent. Sánchez et al. (2013) made efforts to address this concern by observing a sample of 150 local (z < 0.03) galaxies using an integral field spectrograph PMAS/PPAK mounted on Calar Alto 3.5 m as part of the CALIFA survey. In addition to the measurements of resolved H ii regions, Sánchez et al. (2013) determine global estimates for galaxy metallicity in a physically motivated way (at one effective radius), as well as the total SFRs based on Hα. Such a sample, even though relatively small but being free from aperture effects, could prove essential in either strengthening or weakening the status of the M_*–Z–SFR relation. The analysis performed on global measures by Sánchez et al. (2013) concluded that no dependence of MZR on SFR existed. Furthermore, they tentatively explained the apparent presence of this dependence in SDSS data to be due to the aperture affects (their Appendix).

Here we reanalyze Sánchez et al. (2013) data using our preferred methodology: metallicity versus relative SSFR in individual mass bins. The results are shown in Figure 10. Except in the lowest mass bin (panel (A)), the anti-correlation between metallicity (determined by Sánchez et al. 2013 using Pettini & Pagel 2004 calibration of O3N2 method) and SSFR (i.e, the offset from the SF sequence as derived with Sánchez et al. 2013 data) is convincingly present (purple lines show linear unweighted fits). In mass bins centered at log M_* = 10.0 and 10.5, there is only a 5% and 3% probability that the anti-correlation is due to chance (obtained using bootstrap resampling; similar results, 3% and 0.1%, are obtained when measurements are perturbed within the error bars). The positive correlation in the lowest mass bin is not statistically significant (there is a 37% probability that it is due to chance), but it is incompatible with very strong anti-correlation expected at those masses. We point out that Sánchez et al. (2013) sample is incomplete at those masses, and the apparent lack of anti-correlation could potentially be due to the apparent size selection present in CALIFA data set (Walcher et al. 2014). We believe that the reason why this dependence was not detected in the analysis of Sánchez et al. (2013) was because the sample was not split by stellar mass (their Figure 4, bottom right panel). Furthermore, in a small sample that lacks extreme star formers, the trend in metallicities will be relatively modest and therefore difficult to spot on a mass–metallicity plot color-coded by SFR (their Figure 4, lower left panel).

Zoom In Zoom Out Reset image size

We conclude that it is very encouraging that the measurements that avoid the issues of SDSS fibers confirm the MZR dependence on SSFR, and in a degree that is comparable to that using more extensive SDSS data.