A HIGHLY CONSISTENT FRAMEWORK FOR THE EVOLUTION OF THE STAR-FORMING "MAIN SEQUENCE" FROM z ∼ 0–6

At present, several different stellar IMFs are used to derive MS properties. These are usually presumed to be universal—i.e., unchanging with respect to time, current and/or past SFH, metallicity, etc. Current evidence is conflicting: Bastian et al. (2010) claim that the IMF is likely universal, while Kroupa et al. (2013) argue that the IMF becomes more top-heavy (i.e., forming higher fractions of more massive stars) with increasing SFRs. Possible ramifications of this for MS evolution are discussed in Davé (2008), but at present the issue remains unresolved. In this work, we assume a universal IMF. Evolution in the IMF as a function of the SFR could change the derived MS slope, and evolution as a function of redshift could affect our evolutionary fits.

The most common of these IMFs are those of Salpeter (1955), Chabrier (2003), and Kroupa (2001), most commonly (but not universally) integrated from 0.1–100 M_☉. These will be referred to as Salpeter, Chabrier, and Kroupa IMFs, respectively. The assumed IMF impacts both the derived masses and SFRs, leading to variations of up to ∼40% (KE12). At present, there are several different factors used to convert between these different IMFs (E07; S07; C09; K11; Z12; Papovich et al. 2011). We choose the IMF offsets taken from Z12 (also seen in S07 and E07) because they have been calculated recently and assume the same SPS model (BC03) that we standardize to here. These take the form

$$\begin{equation} M_{*,K} = 1.06\,M_{*,C} = 0.62\,M_{*,S}, \end{equation} \tag{ 2 }$$

with the subscripts referring to Kroupa, Chabrier, and Salpeter IMFs, respectively. These correspond to mass offsets of $C_{M_*,C} = +0.03$ and $C_{M_*,S} = -0.21$ dex. These agree well with the SFR offsets used to convert from Kennicutt (1998a; K98) to KE12 (which assume Salpeter and Kroupa IMFs, respectively) for SFRs derived from the FUV and NUV. In all cases, the shift between a Chabrier and Kroupa IMF is essentially negligible. Although all adjustments have been applied for completeness, we note that our results are unchanged if the Chabrier IMF-derived masses are left as they are.

This mass adjustment is functionally equivalent to shifting the MS left or right (i.e., increasing/decreasing the SFR at a given mass) with the observed mass ranges adjusted accordingly (see Tables 4 and 6). These lead to calibration offsets in the normalization, C_M, of $-\alpha \times C_{M_*}$. For a slope of unity, these changes merely result in a shift of the observed range of the MS relation rather than the actual MS relation itself. However, as the majority of data compiled here have slopes of less than unity (see Figure 1), and the SFR offsets are not equivalent to the mass offsets in some cases (see KE12), the majority of these changes do impact the observed normalizations significantly. For these reasons, we choose to only apply explicit IMF adjustments to the derived masses, as the L–ψ relations outlined in the next section implicitly include such adjustments.

SFRs are calculated based on observed galaxy luminosities over spectral ranges that correlate with active star formation in the past 10–100 Myr. These most commonly are the UV continuum (from ∼1500–2800 Å), Hα emission, and the total IR (TIR) continuum (from ∼3–1100 μm). In addition, other SFR indicators, such as 1.4 GHz emission, have further been developed by exploiting the tight observed radio–IR correlation (Condon 1992; Yun et al. 2001; Bell 2003), as well as from SED fitting to individual bands (see S12) or multiband photometry (see M13). These indicators are sensitive to the SFR on different timescales: while Hα probes SFRs on <10 Myr timescales, UV and TIR (and by extension 1.4 GHz) probe SFRs on ∼100 Myr timescales.¹¹ For additional discussion on the nature of SFR indicators and the assumptions used to derive them, see Hao et al. (2011), Murphy et al. (2011, 2012), and KE12.

Most notably, the studies included here calculate integrated luminosities over the entire wavelength range of interest by fitting specific templates to observed bands. In the IR, these templates most often are taken from Chary & Elbaz (2001; CE01), Draine & Li (2001; DL01), Dale & Helou (2002; DH02), and Draine & Li (2007; DL07). In the UV, the most commonly used templates are taken from BC03, although Brammer et al. (2008; B08) and Brammer et al. (2011; B11)¹² are also used. To account for additional strong emission lines, Charlot & Longhetti (2001; CL01) templates are also sometimes used.

Each of these SFR indicators traces ψ in different ways, over different timescales, and with different calibration issues (see, e.g., Table 1 and Section 3 of KE12), with different L–ψ conversions differing by up to ∼50% (see, e.g., the radio SFR calibrations from Yun et al. 2001 and Bell 2003). The standard calibration for most calculated SFRs today is K98, based on a single power-law Salpeter IMF. While K98 gave reasonable SFR calibrations between SFR indicators, for many other wavelengths often studied, the relative calibrations are sensitive to the precise form of the IMF.

KE12 have taken advantage of major improvements in stellar evolution and atmospheric models over the last decade to update the L–ψ relations presented in K98 to a Kroupa IMF, a broken 2-part IMF with a turnover below ∼1 M_☉. A Chabrier IMF, which has a log normal distribution from 0.1–1 M_☉, yields nearly identical results to those of KE12 (Chomiuk & Povich 2011). As KE12 provide a self-consistent set of L–ψ relations for a more realistic IMF (see their Table 1), we opt to convert all previously derived SFRs to this new metric. The ratio of the L–ψ relationships used in individual papers relative to those of K98 and KE12 are listed in Table 3. For SFRs derived from a combination of IR and UV data, we weigh ψ_UV and ψ_IR according to the calibration presented in Section 3.1.5. For more information on these conversions, see KE12, Murphy et al. (2011), and Hao et al. (2011). See Ranalli et al. (2003), Rieke et al. (2009), and Calzetti et al. (2010) for L–ψ conversions at 2–10 keV, 24 μm, and 70 μm, respectively, and Murphy et al. (2012) for an empirical comparison of the radio SFR calibration presented here. Additional composite L–ψ relations (i.e., multi-wavelength dust corrections) can be found in Kennicutt et al. (2009) and Hao et al. (2011). See Calzetti et al. (2007, 2010) for more discussion on many of the issues presented here. For convenience, we include a short description of the SFR calibrations used in this work below.

Assuming a solar metallicity and a constant SFR, Murphy et al. (2011) find that Starburst99 stellar population models yield a relation between the SFR and the production rate of ionizing photons, Q(H⁰), of

$$\begin{equation} \log \psi = \log Q(H^{0}) - 53.14, \end{equation} \tag{ 3 }$$

for Q(H⁰) measured in s⁻¹ and a starburst age of ∼100 Myr. Assuming Case B recombination and an electron temperature T_e = 10⁴ K, the Hα recombination line strength is then related to the SFR via

$$\begin{equation} \log \psi _{{\rm H}\alpha } = \log L_{{\rm H}\alpha } - 41.27, \end{equation} \tag{ 4 }$$

for L_Hα measured in erg s⁻¹. This is a factor of 0.68 that of the corresponding calibration from K98 and probes (0-3-10) Myr (min-mean-90%) timescales. Note that the two coefficients are nearly independent of starburst age for ages ≳ 10 Myr.

As the integrated UV spectrum is dominated by young stars (K98; S07; Calzetti et al. 2005), it is a sensitive probe of recent star formation activity. By convolving the output Starburst99 spectrum with the Galaxy Evolution Explorer (GALEX; Martin et al. 2005) FUV transmission curve, Murphy et al. (2011) find

$$\begin{equation} \log \psi _{{\rm FUV}} = \log L_{{\rm FUV}} - 43.35, \end{equation} \tag{ 5 }$$

for L_FUV measured in erg s⁻¹. This is a factor of 0.63 that of the corresponding calibration from K98 and probes (0-10-100) Myr timescales. Likewise, for the NUV, they find

$$\begin{equation} \log \psi _{{\rm NUV}} = \log L_{{\rm NUV}} - 43.17, \end{equation} \tag{ 6 }$$

for L_NUV measured in erg s⁻¹. This is a factor of 0.64 that of the corresponding calibration from K98 and probes (0-10-200) Myr timescales.

Due to the presence of dust, much of the light emitted by young stars in the UV is absorbed and re-emitted in the IR. In order to derive a calibration for the TIR, Murphy et al. (2011) assume that the entire Balmer continuum is absorbed and re-radiated by dust and that the dust emission is optically thin. After integrating the output Starburst99 spectrum from 912–3646 Å, they find

$$\begin{equation} \log \psi _{{\rm TIR}} = \log L_{{\rm TIR}} - 43.41, \end{equation} \tag{ 7 }$$

for L_TIR measured in erg s⁻¹. This is a factor of 0.86 that of the corresponding calibration from K98 and probes (0-5-100) Myr timescales. Note that the exact timescales are sensitive to SFH (see, e.g., Hayward et al. 2014).

To derive radio SFRs, most studies use the tight, empirical IR–radio correlation (de Jong et al. 1985; Helou et al. 1985; Yun et al. 2001; Bell 2003). This relation is most often expressed in terms of q_IR, where

$$\begin{equation} q_{{\rm IR}} \equiv \log \left(\frac{L_{{\rm IR}}}{3.75 \times 10^{12}\,L_{1.4}}\right), \end{equation} \tag{ 8 }$$

and

$$\begin{equation} L_{1.4} = 9.52 \times 10^{18}\,S_{1.4}\,d_L^2\,4\pi (1 + z)^{-0.2}, \end{equation} \tag{ 9 }$$

where d_L is the luminosity distance of the galaxy in Mpc, S_1.4 is the 1.4 GHz flux density in Jy, and a radio spectral index ($S_\nu \propto \nu ^{\alpha _S}$) of α_S = −0.8 is assumed (e.g., D09; K11). For L_IR ≡ L_TIR, q_IR = 2.64 ± 0.26 dex for SFGs in the local universe (Bell 2003); for L_IR ≡ L_FIR, q_IR is instead 2.34 ± 0.26 dex (Yun et al. 2001). Using the Bell (2003) q_IR value, Murphy et al. (2011) find

$$\begin{equation} \log \psi _{1.4} = \log L_{{\rm 1.4}} - 28.20, \end{equation} \tag{ 10 }$$

for L_1.4 measured in erg s⁻¹ Hz⁻¹. This probes star formation activity in the last ∼100 Myr.

In order to derive masses and SFRs, studies need to assume a specific SPS model. The basic ingredients needed to generate an SPS model are relatively straightforward, and are discussed extensively in Conroy (2013). Unfortunately, systematic uncertainties in calculating particular phases of stellar evolution, inadequacies in current stellar libraries, and other simplifying assumptions can lead to significant errors that are frequently not taken into account (Maraston 2005; Conroy et al. 2009; Percival & Salaris 2009; Behroozi et al. 2010; Conroy & Gunn 2010; Conroy 2013). For example, uncertainties in modeling the little-understood evolution of thermally pulsating asymptotic giant branch (TP-AGB) stars, blue stragglers (BS), and horizontal branch (HB) stars, all of which are relatively luminous, are significant and can have major impacts on the integrated stellar spectrum (Maraston 2005; Melbourne et al. 2012) ranging from ∼0.1–0.3 dex depending on SPS model (Salimbeni et al. 2009; Conroy et al. 2009; M10) in a way that is likely mass-dependent (Salimbeni et al. 2009). SPS calculations also implicitly assume a well-sampled (i.e., fully populated) and unchanging IMF, which may not always be satisfied (Kroupa et al. 2013).

Multiple SPS models are used when fitting for masses and deriving photometric redshifts (photo-zs; see Section 3.2.2). The models used in the compilation presented here¹³ are taken from Fioc & Rocca-Volmerange (1997, 1999) (PEGASE.2), BC03, Maraston (2005; M05), CB07, Charlot & Bruzual (2011; CB11),¹⁴ Polletta et al. (2007; P07), Rowan-Robinson et al. (2008; R08), and Gruppioni et al. (2010; G10). In this study, all masses are calibrated as best possible to BC03 models, as described in Appendix A. PEGASE.2 models are assumed to be similar to BC03 since they use similar stellar evolution tracks (i.e., the Padova 1994 stellar evolution tracks),¹⁵ and so their derived masses are left unchanged. R08 models, although empirically grounded, are regenerated to higher resolution (and given physical parameters) based upon the SPS models Poggianti et al. (2001). These again use similar stellar evolutionary tracks as BC03 models, and so are assumed to be similar. The models of P07 and G10 are fit only in addition to BC03 models in the studies listed here. As the relative rate of their fitting procedure relative to their BC03 counterparts is not detailed in any of the studies provided, possible differences are not accounted for here. Our assumption that these models lead to broadly similar physical parameters (at least for masses) are also supported by M13 at low redshift, who find that using several different SPS models (e.g., BC03, PEGASE) for the same set of priors results in almost identical stellar mass functions.

M05, CB07, and CB11 models, however, utilize different prescriptions to treat the TP-AGB phase that substantially differ from BC03 models. As the TP-AGB phase tends to dominate much of the starlight at certain wavelengths, the revised prescriptions tend to revise masses downward. We treat these models as identical because they implement similar TP-AGB prescriptions (R12), and implement an adjustment upward of $C_{M_*,S} = +0.15$ dex here based on the results of M10 (also an approximate average between the results of Salimbeni et al. (2009) and Conroy et al. (2009)). Most of the adjustments implemented in this way have the fortunate coincidence of being at similar redshifts (z ≳ 2) and being selected via Lyman-break criteria (see Section 3.3.2). The exception is So14, for which the offset is closer to ∼1.6 ($C_{M_*,S} = +0.20$ dex; D. Sobral, private communication). This gives SPS calibration offsets of $C_S = -\alpha \times C_{M_*,S}$.

Some studies choose to eschew using SPS models and SED fitting altogether in favor of analytical M_*/L relations (calibrated on SPS models; e.g., McCracken et al. 2010 and González et al. 2011) which can applied to a wider selection of data to get "cheap" masses (as in P09 and B12). In principle, since these relationships are derived from given SPS models, they should yield good masses on average for similar samples. In addition, many of these M_*/L relations include a built-in color dependence that accounts for variation across the population (i.e., SF versus quiescent; Bell et al. 2003). Extending these relationships to larger samples and a wide range of masses, however, might lead to systematic effects in the derived masses relative to those derived directly through SED fitting. For instance, Ilbert et al. (2010) find that using the analytical M_*/L relationships of Arnouts et al. (2007) for galaxies in the COSMOS (Scoville et al. 2007) field overpredict masses by an average of 0.2–0.4 dex at fixed luminosity. Based on these findings, we adjust P09's masses by an additional −0.2 dex (as they are derived from COSMOS field galaxies, albeit using a slightly different K-band M_*/L conversion), but not those of B12 (which have not been investigated in a similar fashion and also are applied to similar data at similar redshifts).

Different SFHs are needed as inputs to generate the SEDs used to derived galaxy physical properties. The most commonly used of these are declining (D) SFHs, taken from an exponentially-decaying burst with ψ(t) = ψ₀ e^−t/τ, where ψ₀ is the SFR at the onset of the burst (and also the scale-factor used in SED fitting procedures), t is the time since the onset of the burst, and τ is SFR e-folding time. D-SFHs are usually modeled with a wide grid of values ranging from tens of Myr to several Gyr (Maraston et al. 2010). Some fitting procedures modify typical D-SFHs by superimposing random starbursts (DRB-SFHs), usually modeled using a tophat function with a constant SFR and a range of intensities and timescales (S07; M13; So14).

Recent studies, however, motivated by the unphysicality of the extremely short ages often derived with D-SFHs—plus the implied functional form of the MS for galaxies undergoing significant mass assembly—have advocated rising SFHs as better functional fits to the MS than D-SFHs. These have taken several forms: that of exponentially rising (R) SFHs, with ψ(t) = ψ₀ e^t/τ (Maraston et al. 2010; Gonzalez et al. 2012); power-law-rising (RP) SFHs, with ψ(t) = ψ₀ t^α (Papovich et al. 2011; but see Smit et al. 2012); and linearly rising (RL) SFHs, with ψ(t) = ψ₀ + (dψ/dt)t (L11). Frequently, constant (C) SFHs are also used as a go-between for the two options, with ψ(t) = ψ₀ (L12). Studies may also include "delayed-τ" (DT) models, with ψ = (A/τ²)te^−t/τ, with A a normalization constant and the rest of the variables defined as above (St14; see also M13). These models allow the construction of D-SFHs (t/τ ≫ 1) and RL-SFHs (t/τ ≪ 1), as well as several that serve as intermediates between the two.

Given all these current different parameterizations of SFHs, however, it seems that the derived masses are largely independent of the chosen SFH assuming reasonable physical constraints on τ (Maraston et al. 2010; R12; So14). The SFRs from R-SFHs/DRB-SFHs in this scenario, by contrast, can be ∼0.3–0.4 dex (i.e., around a factor of 2) higher than those from simple D-SFHs (Maraston et al. 2010; So14). Often, however, τ values do not choose "physical" values, an artifact of the "outshining" problem—i.e., that the youngest stars tend to dominate the SED in the UV, where the SFR is best constrained, and thus provide poor constraints on stellar ages without extensive multiband photometry and possibly more complex and physically motivated SFHs. Instead, fits frequently choose incredibly small, unphysical values of τ that simply provide the best formal fits to the SED. When τ is free to choose these small values, Maraston et al. (2010) find that the R-SFHs tend to derive masses of ∼0.2–0.3 dex less than those of D-SFHs, and SFRs ∼0.5–0.6 greater than those of D-SFHs.

This problem is not fully rectified by using slightly more complex, two component SFHs (S11; So14), and frequently requires ad hoc limits on τ. Indeed, findings from So14 and Behroozi et al. (2013) seem to indicate that the uncertainty in SFHs leads almost all fitted parameters other than the mass to be extremely unreliable. Given that the SFH of a typical MS galaxy, which likely includes both a rising and declining SFH component of varying degrees as a function of observed mass (and hence formation time; see Section 5.2), as well as the small sample sizes of both studies, we do not utilize any possible SFH-based SFR adjustments here.

Most studies have measured extinction/attenuation photometrically¹⁶ by dust using E(B–V), either derived through SED fitting (E_S) or using the IR-to-bolometric luminosity ratio (IRX) via the IRX-β (where β is the UV slope) relation (E_β) of, e.g., M99. In the literature, E_S values tend to be used cautiously because of the degeneracies between age and reddening (and hence the assumed metallicity and SFH), the parameterization of the extinction curve, and the very limited grid space, leading E_β values to be preferred. However, while observational methods such as the IRX–β relation have for a long time been found to give accurate UV-corrected SFRs compared to those derived from UV+IR observations (B12), it exhibits a significant amount (up to an order of magnitude in some cases) of scatter (Boquien et al. 2012). In addition, results from Wuyts et al. (2011a, 2011b), Price et al. (2013) imply that simple extinction corrections are insufficient to accurately correct for dust, and that more complex geometrical (i.e., patchy) dust models are needed.

Regardless, by observing UV versus UV+IR emission from an ensemble of SFGs, one can apply average extinction corrections that should be sufficient to convert from the observed UV-derived SFR to the bolometric SFR. As observations imply that average IRX of galaxies evolves strongly at higher redshifts (e.g., B12), we will approximate the average IRX observations using results presented by R12a (ψ_bol ∼ (5.2 ± 0.6) ψ_UV) at low-z (z < 4), and B12 (ψ_bol ∼ (2.5 ± 0.5) ψ_UV at z ∼ 4–5; see their Table 6 for the full list of corrections) at high-z (z > 4). We apply these IRX values to weight the corresponding ψ_UV and ψ_IR components from ψ_{UV + IR} data accordingly when adjusting SFR values using the L–ψ relations from KE12 as well as to correct for dust attenuation in data which only reports observed UV luminosities.

Multiple extinction curves have been used in the literature to account for the effects of dust on the observed SEDs in SPS-generated spectra. The ones used in the papers presented here are taken from: Prevot et al. (1984; P84), from observations of the Small Magellanic Cloud (SMC); Cardelli et al. (1989; C89), from various sources in the optical and NIR; Calzetti et al. (2000; C00), from observations of SFGs; Madau (1995; M95) and Charlot & Fall (2000; CF00), from observations of nebular attenuation; and CE01, from observations of local galaxies. A hybrid C00 model with a bump at 2175 Å (C00b) to account for graphite and polycyclic aromatic hydrocarbon (PAH) features is also used in some cases (Ilbert et al. 2009). Although the impact of extinction can be as high as a factor of ∼5 (R12a), the impact of using different extinction curves appears negligible¹⁷ (Papovich et al. 2001; Dickinson et al. 2003). This will not be accounted for here.

However, while the use of different extinction curves might produce similar results, using uniform extinction curves might still produce subtle biases in derived physical results. In particular, if the general shape/amount of PAH emission is correlated with the fitted amount of extinction, then the use of models with constant (or a lack of) 2175 Å features will produce notable biases in SED-derived galaxy properties. Using a flexible parameterization of dust attenuation, Kriek & Conroy (2013) report a negative correlation between the slope of the attenuation curve and the strength of the 2175 Å bump (i.e., SED types with steeper attenuation curves have stronger bumps.). They find this leads to biases in derived dust attenuation (large) as well as masses (small) and specific SFRs (sSFRs, ϕ ≡ ψ/M_*; also small). In addition, they find edge-on and/or low-sSFR galaxies tend to have steeper attenuation curves, while face-on and/or high sSFR galaxies tend to have shallower attenuation curves, implying possible dependencies on orientation. Taking these findings into account may better improve future SED-fitting procedures.

While these findings imply that current SED-fitted physical parameters might display parameter-dependent systematic biases (but see Garn & Best 2010), we do not attempt to account for this effect here.

Strong emission lines, such as Lyα, Hα, Hβ, [O ii], and [O iii] can significantly alter the SED by contaminating observed band photometry. These lines decrease (i.e., make more luminous) the observed magnitude in a given band by up to several tenths of a mag at higher redshifts (Ilbert et al. 2009; S11; Stark et al. 2013). These differences, not accounted for (correctly) by most SPS models (Ilbert et al. 2009), can significantly affect the derived physical parameters taken from the SED fitting process and impact the quality of derived photo-zs. The relative impact depends on the number of bands included in the fit, their respective width, and the redshift of the source: for surveys with a large number of bands (e.g., COSMOS), this effect will be somewhat washed out; however, for surveys with only a handful of bands, this effect can make a big difference (Kriek & Conroy 2013).

Stark et al. (2013) show the effects that emission line contributions can have on the observed M_{UV, 1500}–M_* (i.e., M_*–ψ) relationship at high redshifts (where emission line contamination is most severe), and demonstrate that while on average the slope of the relation remains the same, the overall fitted masses decreases substantially. We choose to implement high-z corrections from their Figure 7 (taken from Robertson et al. 2013), which lead to mass corrections of $C_{M_*,E} \sim (0,-0.03,-0.18,-0.40)$ dex for galaxies at z ∼ (4, 5, 6, 7), respectively. Like Robertson et al. (2013), we have chosen to apply the correction without the hypothesized redshift evolution of the Hα equivalent width (EW) due to the age dependence it would introduce. If we had taken these into account, the corrections listed above would be even larger (e.g., up to an order of magnitude at z ∼ 7). This leads to MS calibration offsets of $C_E = -\alpha \times C_{M_*,E}$.

The effects of differing cosmologies are accounted for by calculating the ratios between luminosity distance, d_L(z), derived from two different cosmologies, and, given the observed redshift range of a sample, applying a $d_L^2$ correction ($C_{d_L}$) at the expected median z of galaxies in the sample after weighting for first-order volume effects. This volume-weighting assumes an approximately constant number density and slowly changing mass distribution of MS galaxies within the redshift bin in question. The effect is negligible in these cases, only changing the derived d_L corrections by less than a percent, but are more significant when we use it later to deconvolve the scatter about the MS (see Section 4.1). Relative to the possible impacts listed in Section 3.1, we find this effect is small, in all cases <0.05 dex. As they are straightforward to derive, however, we choose to apply them out of completeness (see Table 5). Because they boost the luminosity of the entire spectrum, cosmology differences should lead to both increased masses and SFRs. Our calibration offset is then $C_C = (1-\alpha) \times C_{d_L}$.

For the majority of galaxies used in the studies included here, redshifts have been derived photometrically (photo-zs) via SED fitting rather than spectroscopically (spec-zs). SPS models are used to derive these photo-zs, which simultaneously provide the masses (and sometimes SFRs) used in these studies. Photo-zs have varying precision, ranging from 0.8%–3% scatter compared to their spec-z counterparts (Ilbert et al. 2013), and can be subject to "catastrophic failures" where the photo-zs and spec-zs disagree by more than 15% (η ≡ |z_phot − z_spec|/(1 + z_spec) > 0.15). Note that these statistics are only available when spec-zs are available, and thus are often based on only the brightest galaxies (which are often targeted in I-band selected surveys).

Besides just misfits caused by bad photometry, a small number of bands, or a multi-peaked redshift probability distribution function (PDF), catastrophic errors can occur systematically by, e.g., confusing the Lyman break at ∼1220 Å and the Balmer/4000 Å break (St14). Although the errors from the average scatter are small, the effect of catastrophic failures on the M_*–ψ relation relative to that of confirmed spectroscopic samples has yet to be fully investigated. We conduct a simple experiment to qualitatively assess the effects of catastrophic errors M_*–ψ relationship, and find that their effect on the overall distribution appears small, even for a large fraction of catastrophic errors (see Appendix B).

In many cases, photo-zs have not been compared with spec-zs across the full mass and redshift ranges to which they have been applied; this serves to both check their accuracy and are often necessary for calibration purposes (Hildebrandt et al. 2010; Abdalla et al. 2011; Dahlen et al. 2013). Existing spec-z or narrow-band selected studies, however, seem to agree well with photo-z derived distributions (e.g., S12; C14). In addition, simulated errors and catastrophic failure rates agree with the measured spectroscopic samples and are accounted for in some works (e.g., Ilbert et al. 2013). Finally, the quality of photo-zs have been checked with pair statistics and cross correlations, which seems to confirm errors derived from spec-zs (Benjamin et al. 2010).

On the whole, photo-z methods do not seem to display large redshift biases relative to spec-zs, and the likely induced scatter is small relative to scatter about the MS and other systematic errors (Hildebrandt et al. 2010; Abdalla et al. 2011; Dahlen et al. 2013). Any possible systematic offsets they have relative to spec-zs are not accounted for in this study.

Besides the variations in generating SEDs that have been detailed above, the SED fitting procedure used to derive photometric redshifts differs for different codes.¹⁸ Each of these fitting procedures, besides contamination from catastrophic errors, might exhibit biases in the determined photo-zs relative to the true spec-zs and/or each other. In particular, the best-fit parameters derived from SED fitting tend to be sensitive to small changes in parameter space and errors on the photometry. This can be reduced by incorporating a wider range of parameter space into the final mass, such as by taking the median mass across all solutions in the entire multi-dimensional parameter space for each fit that lies within 1σ of the best fit (So14; see Appendix E). At the moment, however, since such procedures are not widely used, we will not attempt to account for these effects here.

In order to directly test different photo-z codes/fitting procedures against each other, Hildebrandt et al. (2010), Abdalla et al. (2011), and Dahlen et al. (2013) compare photo-z code performance against each other using identical samples. Their results indicate that, in general, all codes produce reasonable photo-z estimates in both an absolute and relative sense, although using a training set of spec-z priors reduces both the scatter and the fraction of catastrophic errors. Their findings also indicate that using a training set from a small region of the sky does not seem to produce biases when applied to larger survey areas, and that the median of all codes seems to do better than any individual code at matching spec-zs.

Most crucially, Dahlen et al. (2013) find that photo-z errors and the fraction of catastrophic errors are the largest for data at higher magnitudes (i.e., are fainter with larger error bars), which implies the majority of photo-z errors should happen preferentially to low-mass, low-SFR galaxies observed within any given sample (precisely where spec-zs are lacking). While these results provide areas for photo-z codes to improves and that should be investigated, based on these overall positive results, we do not opt to attempt to account for possible differences among SED fitting procedures.

Stellar evolutionary tracks (i.e., isochrones) used by SPS models can be strong functions of metallicity. Currently, SPS models do not model metallicity evolution self-consistently, which would involve tracing the evolving metallicity content of stellar populations over time from supernovae injections, mixing, elemental abundance patterns, etc., and their subsequent impact of star formation and evolution. Instead, many resort to using simple stellar populations (SSPs), which follow the evolution in time of the SED of a single, coeval stellar population at a single fixed metallicity and abundance pattern. The effects of using SSPs relative to populations where metallicity evolution is taken into account is not fully understood. As SSPs are utilized in all SPS models considered here and a fundamental assumption in the derivation of physical parameters from SEDs, we take this to be an unknown systematic that cannot be quantified and/or accounted for at this time. See Conroy (2013) for further discussion.

While systematics from using fixed-metallicity SSPs are not accounted for, we can at least investigate a related assumption: the effects using different metallicities in SSPs have on the derived physical parameters. At low redshifts (z ≲ 1), results from M13 (see their Appendix B) seem to indicate that assuming a fixed solar metallicity relative to a much wider metallicity distribution (from ∼0.2–1.5 Z_☉) does not have a major impact on the resulting mass distribution both at fixed redshift and as a function of redshift (their impact on SFRs has yet to be thoroughly investigated). Based on these findings, and the fact that the majority of studies included here include sensible metallicity priors, we do not attempt to implement any adjustments due to possible metallicity-induced effects.

While the mean and median of a log-normal distribution are approximately identical when calculated in log space, the expected mean of a log-normal distribution is skewed in linear space (Behroozi et al. 2013). This offset depends on the scatter present in the distribution—for a log-normal distribution with a median of 1 and scatter σ (dex), the expected mean will instead be

$$\begin{equation} \langle x \rangle = \exp [0.5(\sigma \ln 10)^2]. \end{equation} \tag{ 11 }$$

This leads to an offset between the mean and median in log space of

$$\begin{equation} \Delta x \approx 1.15\sigma ^2. \end{equation} \tag{ 12 }$$

For an intrinsic scatter of σ = (0.30, 0.35) dex, this corresponds to Δx = (0.10, 0.14) dex. As all radio data included here (D09; P09; K11) have used median stacks to find the true mean of the SFR for a given mass bin (White et al. 2007), this effect translates to a systematic overestimation of the SFR at a given mass by approximately 0.1 dex compared to most other data included here. This effect has been included in the C_ψ calibration offsets presented in Table 5.

Selection effects within each study—not to mention within the definition of the MS itself—also can affect both the derived slopes and their evolution as a function of redshift (O10; K11; W12). While most of these (see Table 3) are efficient at selecting SFGs, they do not all select the same population. As K11 show in their Appendix C, B–z versus z–K (sBzK; Daddi et al. 2004a (D04)) selection—and a bluer selection criteria in general—is biased toward more "active" SFGs (i.e., with higher (s)SFRs), excluding good portions of galaxies that are classified as SFGs via other selection mechanisms (NUV − r versus r−J; NUVrJ), and give steeper MS slopes (see also O10 and K11). W12 shows that this effect further translates into an inherent bias against redder, more dust-attenuated SFGs (see their Figure 3), which have lower slopes compared to their bluer, less dust-attenuated counterparts. See also Sobral et al. (2011) for more discussion on this issue.

Because of these effects, selection methods that are inherently biased toward bluer, highly active, non-dusty SFG populations will preferentially select a subset of the MS population with a higher slope relative to other selection mechanisms. These selection methods include: sBzK, used to select SFGs from 1.4 < z < 2.5 (D07; R11; K13); the Lyman break (Steidel et al. 1999; Stark et al. 2009; Bouwens et al. 2011; B12), used to select high-z Lyman-break galaxies (LBGs); and U − g versus M_bol, or any other cut on the color–magnitude diagram (CMD) that explicitly selects based on (blue) color (E07). We therefore classify these methods as "bluer" selection mechanisms, along with luminous infrared galaxy (LIRG) selected samples such as that of E11—although LIRGs are definitely SFGs, nearby LIRGS tend to be highly active SFGs with large amounts of dust attenuation and extreme amounts of star formation, in contrast to "regular" MS galaxies that are more similar to the Milky Way at low redshifts from, e.g., B04 and S07.

As can be seen by comparing selection methods from, e.g., E07 (see their Figure 2) and Ilbert et al. (2013) (see their Figure 3), non-"bluer" selection methods (broadly classified as "mixed") seem to provide not only a "cleaner" cut between SFGs and quiescent galaxies, but a more diverse star-forming population. The total classification scheme of these two selection types is listed in Table 3. As mentioned earlier, while these different selection methods do not seem to affect the average observed SFRs across different publications, they do seem to influence the derived slopes and the intrinsic scatter (see Section 3.3.3).

Although differences between bluer and mixed selection criteria can lead to differences in the derived MS relations, all MS studies should ideally only include SFGs in their analysis. Several of the studies included here do not opt to impose a color–color cut of some sort to separate out SFG and quiescent galaxy populations¹⁹ (C09, So14, and C14). These "non-selective" studies consequently display prominent differences from SFG-only studies. As quiescent galaxies "contaminate" their highest mass bins at a wide range of redshifts, their lower SFRs significantly reduce the slope. In addition, their increased prevalence at lower masses at lower redshifts (as more and more galaxies "quench") leads to increasing offsets in normalizations for any flux-limited survey. This effect is accentuated by increases in survey sensitivity, which can drive the SFR floor lower at all masses.

As expected, we find that all non-selective studies agree with other data relatively well at lower masses (especially at higher redshifts), but disagree significantly at higher masses due to shallower MS relations.²⁰

In Figure 1, we plot the derived slopes of each MS sample color-coded by selection method of the parent sample, the subsample used for analysis, and our groupings listed above. We find that, while some biases in MS selection might emerge from parent samples selected primarily on restframe UV, the majority of biases occur in the precise selection of the subsample. As expected, we also find that "bluer" MS observations display slopes between ∼0.8–1 and are relatively similar over the majority of the age of the universe, while "mixed" observations center around ∼0.6 and display possible time-dependencies.²¹ Based on these results, we decide to use our bluer/mixed classification scheme to account for different biases inherent in SFG/MS selection, preferring mixed selection methods to bluer ones since they give us a larger and more diverse SFG sample while still excluding most quiescent galaxies.

We find that the true and deconvolved scatters (see Section 4.1) reported in all sBzK-selected studies (D07; R11; K13) are systematically lower than reported in other papers, even for large sample sizes (R11). Furthermore, their resulting values are low enough to likely be unphysical, especially given the possible ∼0.1 dex of intrinsic scatter in determining mass (even when considering possible convariances; see Appendix E). This seems to indicate that the scatter observed in these papers is not representative of the redshift range that they encompass, and hence that sBzK is substantially biased compared to other selection mechanisms (data taken from LBGs, for instance, show similar scatters as other "mixed" samples; M10; L12; R12).

Most likely, this difference is due to an inherent bias built into the sBzK selection mechanism itself. As outlined in D04, the B–z/z–K line used to select sBzK galaxies was designed to be parallel to the reddening vector. However, due to the age-extinction degeneracy, this means that age runs perpendicular to the selection function, and implies that you will systematically be missing older (and hence likely more massive and dusty) galaxies because of their redder colors. As LBG selection mechanisms are not explicitly designed this way, although redder SFGs are still selected against, the selection bias is not as systematic or complete as using sBzK. Thus, while sBzK is effective at selecting for SFGs for 1.4 < z < 2.5, the distribution and scatter of the sample is biased (σ_{t, BzK} ∼ 0.1 while σ_{t, med} ∼ 0.2; see Section 4.1.2).

However, tests we have conducted on the Ilbert et al. (2009) COSMOS catalog find that sBzK-selection is actually quite efficient at selecting out SFGs between 1.4 < z < 2.5 (O. Ilbert 2013, private communication; although see D09, P09, and K11). While this does not rule out the possibility of intrinsically biased selection, it does point to another possibility. Instead, the narrower distribution might likely arise due to systematic biases in calculating bolometric SFRs. Extinction corrections in sBzK samples are determined exclusively from B–z color (D04), and so the same bands used to select the sample are also used to determine the dust attenuation. Using a large sample of LBGs, B12 finds that using the same passbands for both selection and dust attenuation measurements leads to large biases in the derived dust attenuations. This implies we might be witnessing a similar problem with sBzK-selected sample here, where the problem is not with inherent selection biases, but with substantially biased calculations of dust attenuation and hence a narrower distribution bolometric SFRs.

We note that apart from sBzK, the only other data points which display lower-than-average scatter are those with small sample sizes (N ≲ 250), as well as those from S09 (although S09 uses a 2σ-clipped fitting procedures that biases the scatters by default; see Table 4).

There are several main biases that characterize observations of the SFG MS: incompleteness, Malmquist bias, and Eddington bias. If a survey is not mass-complete, observations will be biased toward bluer galaxies both due to the flux-limited nature of most surveys as well as many SED fitting procedures (which are sensitive to the signal-to-noise ratio of the photometry; see Dahlen et al. 2013), which will affect properties of the MS relation below the mass completeness limit. This is easily rectified by only using data where the survey is approximately mass complete, as is done in, e.g., W12 (see their Figure 1). For the studies collected here, we find that on average this effect on the reported MS fits is small compared to the other issues discussed above, and therefore do not correct for it here.

There are also other competing effects in most surveys that tend to become more prominent near mass limits, most notably a Malmquist bias of selecting galaxies with larger SFRs at a given stellar mass in a flux-limited sample (see R12's Appendix B). As the (s)SFR of galaxies are a strong function of their mass, Malmquist bias will result in higher (s)SFRs derived on average for a given mass for masses where the flux limit approaches the hypothesized distribution. As shown in R12 (see their Figure 26), this bias can lead to derivations of sSFRs from their true values on the low mass end by a factor of ∼3–4. In general, such a bias is strongest the more flux-limited a sample is (and as such is different from just strict mass completeness), and most prominently affects galaxies on the low-mass end of the MS. As this leads to higher average SFRs for these objects as compared to higher mass objects, this would lead to a shallower fitted slope.²²

Another possible impact of Malmquist bias would be a strong selection effect toward bursting low-mass SFGs near the detection limits, as they would be more likely to be detected over their non-bursting counterparts. This might lead to a strong bias in SED-fitting procedures toward very young ages for low-mass, low-SFR systems (bottom-left of the MS), which might in turn bias MS slopes. As such a trend is not seen in R12, such a strong systematic bias is likely not strong.

On the high-mass/high-SFR end, Eddington bias (i.e., that random scatter in a given mass/luminosity bin will preferentially scatter objects up into higher mass/luminosity bins because they have comparatively fewer objects) tends to be much more dominant. Such bias might lead to a flattening of the MS at high masses as lower mass (and hence lower SFR) objects are scattered up into higher mass bins. This effect, however, should also lead to objects with lower SFRs being upscattered into higher SFR bins (assuming, of course, that the SFR/mass derivations are somewhat independent of one another, as they are constrained by different portions of the SED), causing an upturn in the MS relation. These two effects should then combine to produce a more densely populated upper population in the derived MS relation, with a downturn for samples with well-constrained SFRs and less well-constrained masses (e.g., empirically derived Hα SFRs from high-quality spectra and SED-fitted masses from multiband photometry), an upturn for samples with well-constrained masses and poorly constrained SFRs (i.e., extinction-corrected UV SFRs with SED-fitted masses from extensive, high-quality multiband photometry), and a similar slope for samples with about equivalent constraints on both (i.e., both masses and SFRs derived through SED fitting).

All these scenarios are only relevant, however, assuming photo-z accuracy does not depend on other physical parameters and is relatively good for the majority objects included in the fit. This is not necessarily true—like masses and SFRs, photo-z accuracy is sensitive to the overall shape of the SED. This can lead to complex covariances which have not been fully explored (but see Appendix B). If the photo-z is incorrect, then masses may likely be overestimated and SFRs derived from the incorrect portion of the spectrum, which will lead to more complicated behavior. Note that not all studies are affected by photo-z biases: sBzK-, LBG-, and line-selected samples have precise redshift distributions that are applied to the mass fitting rather than using the photo-z for each individual source.

In all cases, however, there will be a bias toward higher numbers of mass/SFR objects. These effects should not have a large impact on MS relations derived directly from nonstacked data (due to the small number of objects at the high mass end), although for stacked data (especially using mean instead of median stacks), this effect is expected be more prominent. In many studies, however, the slope of the MS is computed from binned data, with mass bins often assigned equal weight regardless of the size of each bin. In these cases, using mass bins is essentially the same as stacking and implies likely contributions from Eddington bias on the high-mass end of the MS. Although we do not attempt to correct for it here, we thus cannot exclude a significant contribution from Eddington bias for massive galaxies.