A TURNOVER IN THE GALAXY MAIN SEQUENCE OF STAR FORMATION AT M_* ∼ 10¹⁰ M_☉ FOR REDSHIFTS z < 1.3

We interpolate the 90% stellar mass completeness limits from Ilbert et al. (2013) to determine approximate mass completeness thresholds at all redshifts, and select only galaxies with stellar masses above their redshift dependent mass completeness limit. When studying star-forming galaxy populations, we separate "star-forming" and "quiescent" galaxies using a modified version of the Williams et al. (2009) two-color selection technique: NUV − r⁺ versus r⁺ − J (as described in Ilbert et al. 2013). Specifically, galaxies with absolute magnitude colors M_NUV − M_r > 3(M_r − M_J) + 1 and M_NUV − M_r > 3.1 are considered "quiescent" (∼15% of the sample) while the remaining galaxies are considered actively star-forming galaxies.

Star-forming galaxies that also contain luminous active galactic nuclei (AGNs) are a concern because the luminosity from the AGN is extremely difficult to separate from emission from star formation, and thus these sources may have erroneously high SFRs (although this concern is lessened for FIR sources because AGNs generally heat dust to temperatures too hot to radiate in the far-infrared). On the other hand, many of the galaxies that host AGNs also contain significant star formation, and removing these sources introduces a bias to our study. We find that the overall results of our study are not significantly affected by either the inclusion or exclusion of these sources, so we do not remove galaxies that have been detected in the X-ray (∼0.5% of sample) by XMM-Newton (Brusa et al. 2010) or Chandra (Civano et al. 2011), or that have IRAC power-law colors (Donley et al. 2008) that suggest AGN activity (∼25% of sample).

The Herschel-selected sample of galaxies is described in detail in Lee et al. (2013, hereafter L13) and is briefly summarized here. L13 use Spitzer 24 μm and VLA 1.4 GHz priors to find 4,218 sources in COSMOS that were each detected in at least two of the five available Herschel PACS (100 μm or 160 μm) and SPIRE (250 μm, 350 μm, or 500 μm) bands. These sources span log(L_IR/L_☉) = 9.4–13.6 and z = 0.02–3.54. Dust properties of each source (e.g., L_IR, T_dust, M_dust) were measured by fitting the full infrared photometry to a coupled modified blackbody plus mid-infrared power law using the prescription given in Casey (2012) and assuming an opacity model where τ = 1 at 200 μm.

There is a population of galaxies that are classified as "quiescent" from their NUV − r⁺ versus r⁺ − J colors, but have been detected in the infrared by Herschel or Spitzer, suggesting that these galaxies are actually undergoing a significant amount of star formation (∼7% of the "quiescent" population). These galaxies are more consistent with being very dusty objects that have extremely red colors due to obscuration, not lack of active star formation, so we include these galaxies in our sample of star-forming galaxies (see Section 4.3).

Ilbert et al. (2013) measure accurate 30 band photometric redshifts of the full K_s-band COSMOS catalog. We find a median Δz/(1 + z) = 0.02 in our sample of star-forming galaxies, with a catastrophic failure (|Δz|/(1 + z) > 0.15) in 5.6% of sources. In addition, physical parameters such as stellar mass and star formation rate have been calculated by fitting the spectral energy distributions (SEDs) to synthetic spectra generated using the stellar population synthesis (SPS) models of Bruzual & Charlot (2003). We also recalculate the physical parameters using different templates and extraction parameters, and find that our final results are not affected by the specific choice of template. Thus, we use the same set of parameters used to create the catalog in Ilbert et al. (2013), but with updated near-IR photometry from SPLASH.

There are many methods for estimating a galaxy's SFR based on observations at various wavelengths (for a review, see Kennicutt 1998; Murphy et al. 2011). Here we compare a few commonly used SFR indicators using a subset of COSMOS galaxies to determine how much the different SFR methods disagree. In all cases, we measure the total SFR as SFR_Tot = SFR_IR + SFR_UV.

3.1.1. Infrared Derived SFR

As discussed in Section 2.2, the infrared properties of the Herschel-selected galaxies have been measured by fitting the infrared SEDs to a coupled modified blackbody plus mid-infrared power-law model (Casey 2012). It has been shown that measuring the L_IR from fitting the far-infrared data to libraries of SED (e.g., Chary & Elbaz 2001; Dale & Helou 2002) gives roughly the same results as the modified blackbody plus power-law model (Casey 2012; U et al. 2012; Lee et al. 2013). The infrared observations give us an estimate of the obscured SFR, and we combine this with UV observations of the unobscured SFR to derive the total SFR as in Arnouts et al. (2013):

$$\begin{equation} {\rm SFR}_{\rm {Total}} = (8.6 \times 10^{-11}) \times (L_{\rm {IR}}+ 2.3 \times \nu L_{\nu }(2300\, \rm {\mathring{\rm{A}}})) \end{equation} \tag{ 1 }$$

where L_IR ≡ L(8–1000μm) and all luminosities are measured in units of L_☉. For sources with SFR ≳ 50 M_☉ yr⁻¹, the infrared contribution dominates the total SFR, contributing as much as ∼90% of the total SFR.

While we have excellent Herschel coverage of the full 2 deg² COSMOS field that yields 4218 sources, the detection limits of Herschel introduce a selection bias against all but the most luminous infrared sources. A common method of determining the L_IR of less luminous galaxies is to use deep Spitzer 24 μm data to estimate the far-infrared luminosity (e.g., Kennicutt et al. 2009; Rieke et al. 2009; Rujopakarn et al. 2013). COSMOS has extremely deep coverage at 24 μm and Le Floc'h et al. (2009) provide SFR estimates for 36,635 galaxies, which we use to extend our sample of infrared detected galaxies to more moderate luminosities.

3.1.2. Optical- and UV-based SFR Indicators

We compare commonly used SFR indicators using a common subset of COSMOS galaxies to determine how much agreement there is between the different measures of SFR. We have a large set of Herschel detected galaxies from Lee et al. (2013) where, for the first time, we have direct measurements of both the obscured and unobscured SFR (from UV observations) at a wide range of redshifts. We compare the other SFR indicators discussed previously (24 μm, SED fits, BzK, and NRK) to this sample of 4218 Herschel detected galaxies.

Figure 1 displays the comparison of SFR_Total from the four different indicators discussed above to SFR_Total as measured by Herschel. Density contours show the location and concentration of the majority of the sources, with outliers shown in gray circles. Median values in 20 equally populated bins of SFR_{Total, Herschel} are overplotted to show average trends. To determine the strength of the correlation between each SFR indicator and SFR_{Total, Herschel}, we measure the Pearson correlation coefficient (ρ) and provide these values at the top of each subpanel. The Pearson correlation coefficient can vary between +1 and −1, with +1 indicating total positive correlation, 0 indicating no correlation, and −1 indicating total negative correlation. We also measure the median difference between each SFR indicator and Herschel SFR (〈Δlog (SFR)〉) and list these values at the top of each subpanel.

Zoom In Zoom Out Reset image size

The 24 μm determined SFR correlates with the Herschel SFR very well (ρ₂₄ = 0.88, 〈Δlog (SFR₂₄)〉 =0.12), except at the highest IR luminosities. This trend has been previously explored in many studies which find that at moderate redshifts and IR luminosities, 24 μm observations are a good proxy for L_IR, but at high redshifts and infrared luminosities, the 24 μm estimates tend to overpredict the true L_IR, possibly due to redshifting of the observed 24 μm band to wavelengths contaminated by polycyclic aromatic hydrocarbon features (e.g., Papovich et al. 2007; Lee et al. 2010; Elbaz et al. 2011). As we are using the 24 μm SFRs to fill in the low and moderate luminosity galaxies that are not detected with Herschel, this discrepancy is not a major issue for our work.

The NRK SFRs also show strong correlation with the Herschel-derived SFRs (ρ_NRK = 0.79, 〈Δlog (SFR_NRK)〉 =0.17). This is not completely unexpected since the NRK method was developed using 24 μm derived SFRs as a baseline, but the NRK measured SFRs match very well with those derived from Herschel. Like with SFR₂₄, the correlation shows signs of breaking down at the highest SFRs, but as long as the NRK is used mainly for low SFR galaxies, it provides a reliable estimate of the SFR.

By contrast, the agreement between SFR_SED and SFR_Total is quite poor, showing much weaker correlation between the two indicators (ρ_SED = 0.56). The tightness of the correlation is also much broader (〈Δlog (SFR_SED)〉 =0.43), even at low SFRs where the median points lie closer to the unity line. Again, at high SFRs the median points show a clear deviation from unity. Wuyts et al. (2011) are able to find a better match between SFR_SED and SFR₂₄ if they tune key parameters of the SED fit, such as τ_min, the e-folding time of the exponentially declining star formation history. The exact tuning needed varies based on several other assumptions in the SED fitting procedure, such as different stellar population synthesis codes, and even when the tuning is done, the computed SFR_SED still systematically underestimates the true SFR for a significant fraction of sources.

Finally, the comparison between SFR_BzK and SFR_Total shows essentially no correlation (ρ_BzK = −0.11). It should be noted that the redshift range of the BzK indicator (1.4 < z < 2.5) limits us to a small sample size containing only the brightest galaxies. As seen in Figure 2, this selection limits our comparison to galaxies at SFRs where all indicators begin to deviate significantly from SFR_Total. Stacking analyses suggest a stronger correlation between average SFR_BzK and average SFR_Herschel at fainter luminosities (Rodighiero et al. 2014), but the tightness of the distribution is not well determined. The BzK galaxies that are Herschel detected show no correlation between the SFR derived from the BzK method and from Herschel measurements.

Zoom In Zoom Out Reset image size

3.2.1. Selection Effects of SFR Indicators

3.2.2. A Ladder of SFR Indicators

While all three non-infrared based SFR indicators fail to accurately estimate the SFR in high luminosity galaxies, the NRK method provides the most accurate and consistent estimates across the full dynamical range of Herschel SFRs. At high SFRs (SFR ≳ 30 M_☉ yr⁻¹), 70% of our sample is directly detected in the infrared by either Herschel or Spitzer. Thus, we can study the full population of star forming galaxies by constructing a "ladder" of SFR indicators (as in Wuyts et al. 2011) based on the Herschel, Spitzer 24 μm, and NRK SFR indicators. All sources have SFR_Total calculated using Equation (1), with different methods of determining L_IR. For sources detected by Herschel, we measure L_IR from fitting the far-infrared photometry to the Casey (2012) graybody plus power-law models. We use the L_IR estimated from 24 μm (Le Floc'h et al. 2009) for sources that are not detected by Herschel but are detected at Spitzer 24 μm. And for the remaining sources, we estimate the L_IR using the NRK-derived IRX (as discussed in Section 3.1.2). Although we include NRK-derived IRX for all galaxies above our mass-completeness limits, the method has only been well-calibrated for M_* > 10^9.3 M_☉. Infrared stacking suggests that any systematic offsets should be small, but when calculating MS relationships we only include galaxies with M_* > 10^9.3 M_☉.

The relative fraction of sources with SFRs measured from each indicator is plotted in Figure 3 as a function of both stellar mass and redshift. Below M_* ≲ 10^9.5 M_☉, SFRs are almost all determined from NRK, but at higher stellar masses, the fraction of sources with direct infrared measurements increases until about 25% (60%) of sources at M ≳ 10^10.5 M_☉ have SFRs determined from Herschel (Spitzer 24 μm). Because of the redshift limitations of the NRK method (see Arnouts et al. 2013) and the larger errors associated with the SFR_SED and SFR_BzK indicators, we restrict the rest of our analysis to redshifts 0 < z < 1.3, where we can more accurately measure the SFR of our full sample.

Zoom In Zoom Out Reset image size

From Figure 4, it is clear that a single power law does not accurately describe the relationship between stellar mass and star formation rate. We split our sample of star-forming galaxies into six equally populated redshift bins, each of which are then split into 30 equally populated stellar mass bins (with ∼350 sources in each bin), and calculate the median SFR in each bin. We limit our sample to stellar masses above a conservative mass limit (see Table 1) to ensure that we are not affected by systematics. The specific number of redshift bins does not affect the following results, although we must balance between having redshift bins that are too wide and combine galaxy samples at different epochs with having redshift bins that are too narrow and are affected by small number statistics. The same is true for the number of stellar mass bins, although having at least 30 bins is preferable for accurately determining the goodness of fit to the models. The median SFRs for every bin are plotted in Figure 5, colored by redshift, with bootstrapped errors on the median represented by vertical bars.

Zoom In Zoom Out Reset image size

Table 1. Best-fit Parameters for SFR/M_* Correlation

Redshift Range	〈z〉	log (Mlimit)	$\mathcal {S}_{0}$	$\mathcal {M}_{0}$	γ	Reduced χ2
		log (M☉)	log (M☉ yr−1)	log (M☉)
0.25–0.46	0.36	8.50	0.80 ± 0.019	10.03 ± 0.042	0.92 ± 0.017	1.74
0.46–0.63	0.55	9.00	0.99 ± 0.015	9.82 ± 0.031	1.13 ± 0.033	1.52
0.63–0.78	0.70	9.00	1.23 ± 0.016	9.93 ± 0.031	1.11 ± 0.025	1.48
0.78–0.93	0.85	9.30	1.35 ± 0.014	9.96 ± 0.025	1.28 ± 0.034	1.84
0.93–1.11	0.99	9.30	1.53 ± 0.017	10.10 ± 0.029	1.26 ± 0.032	0.62
1.11–1.30	1.19	9.30	1.72 ± 0.024	10.31 ± 0.043	1.07 ± 0.028	1.24

Notes. Parameters of the best-fit model to the star-forming main sequence. The full sample of 62,521 star-forming galaxies is split into six equally populated bins, with each bin containing ∼17, 745 galaxies. Within each redshift bin, the galaxies are split into 30 equally populated bins of stellar mass. The median SFR in each mass bin is calculated and then fit to $\mathcal {S} = \mathcal {S}_{0} - \log\, [ 1 + ({10^{\mathcal {M}}}/{10^{\mathcal {M}_{0}}})^{- \gamma } ]$, where $\mathcal {S} = \log ({\rm SFR})$ and $\mathcal {M} = \log (M_{*})$. Table columns are as follows: (1) redshift range of bin; (2) median redshift; (3) stellar mass limit of redshift bin; (4) $\mathcal {S}_{0}$, the maximum value of $\mathcal {S}$; (5) turnover mass; (6) low-mass power-law slope; (7) reduced χ² of fit.

Download table as:  ASCIITypeset image

Using the MPFIT package implemented in IDL (Markwardt 2009), we fit the median log (M_*) and log (SFR) in each redshift bin with many models including linear, second order polynomial, and broken linear, and find the best fit is provided by the following model:

$$\begin{equation} \mathcal {S} = \mathcal {S}_{0} - \log \left[ 1 + \left(\frac{10^{\mathcal {M}}}{10^{\mathcal {M}_{0}}} \right)^{- \gamma } \right] \end{equation} \tag{ 2 }$$

where $\mathcal {S} = \log ({\rm SFR})$ and $\mathcal {M} = \log (M_{*}/M_{\odot })$. We choose this model because (1) at all redshifts, it provides the best reduced χ² fit to the data, and (2) unlike polynomial fits, the parameters of the model allow us to quantify the interesting characteristics of the relation between stellar mass and SFR: γ, the power-law slope at low stellar masses, $\mathcal {M}_{0}$, the turnover mass (in log (M_*/M_☉)), and $\mathcal {S}_{0}$, the maximum value of $\mathcal {S}$ (or the maximum value of log (SFR)) that the function asymptotically approaches at high stellar mass. The best-fit parameters for each redshift bin are listed in Table 1.

In the top panels of Figure 6, we plot the evolution of $\mathcal {S}_{0}$, $\mathcal {M}_{0}$, and γ as functions of log (1 + z). The bottom panels of Figure 6 examine the covariance between these parameters by displaying the 95% confidence error ellipses.

Zoom In Zoom Out Reset image size

We see clear and strong evolution of $\mathcal {S}_{0}$ with redshift, and the best fit line suggests an evolution of $\mathcal {S}_{0} \propto (4.12 \pm 0.10) \times \log (1+z)$, or equivalently, SFR₀∝(1 + z)^{4.12 ± 0.10}. The covariance between $\mathcal {S}_{0}$ and $\mathcal {M}_{0}$ (Figure 6(D)) and between $\mathcal {S}_{0}$ and γ (Figure 6(E)) is relatively minor, so we infer that the evolution in $\mathcal {S}_{0}$ is true evolution and not due to variation in the other parameters.

Both $\mathcal {M}_{0}$ and γ show some evidence for weak evolution to more massive $\mathcal {M}_{0}$ and steeper γ with redshift, although much of the perceived evolution may be due to the covariance seen in Figure 6(F). The best (linear) fit to the evolution in the turnover mass is given by $\mathcal {M}_{0} \propto (1.41 \pm 0.20) \times \log (1+z)$. We test the possible redshift evolution of turnover mass by calculating where the data deviates by 0.2 dex from a single power law fit to the low mass data and find similar evolution, suggesting that the turnover mass does indeed change with cosmic time. The low-mass power-law slope, γ, has a best-fit line that suggests evolution of γ∝(1.17 ± 0.13) × log (1 + z). Redshift evolution in γ to steeper slopes at earlier cosmic times would suggest that the SFR in the lowest stellar mass galaxies does not increase as much as in more massive systems.

We have described the MS relationship between SFR and M_* for star-forming galaxies. However, possible misclassification of galaxies as either "star-forming" or "quiescent" could drastically affect the trends we observe.

As described in Section 4, we remove galaxies that are considered quiescent from our sample using the selection M_NUV − M_r > 3(M_r − M_J) + 1 and M_NUV − M_r > 3.1 (Ilbert et al. 2013). This selection is shown in Figure 7, with the full mass-selected sample of COSMOS galaxies at 0.2  <  z  <  1.3 generally separated into two distinct "star-forming" and "quiescent" regions. Improper classification of the "in between" galaxies that are not obviously star forming or quiescent could lead to changes in the MS shape. To test this possibility, we shift the entire separating line (both horizontal and diagonal segments) between quiescent and star-forming galaxies by ±0.4 mag in M_NUV − M_r, and in either case there is no appreciable change to the MS.

Zoom In Zoom Out Reset image size

Galaxies detected in the infrared by Herschel or Spitzer 24 μm are highlighted in Figure 7, and while the majority fall on the star-forming sequence, we see a number of objects that lie in the quiescent region. This population of infrared-detected quiescent (IR-Q) galaxies is relatively small, with only ∼7% of the galaxies classified as quiescent having detectable infrared emission, but these misclassified galaxies are predominantly found at high stellar mass. The fraction of quiescent galaxies detected in the infrared increases rapidly from ⩽1% at M_* ∼ 10^9.5 M_☉ to 15%–20% at M_* ⩾ 10¹¹ M_☉, and this trend holds at all redshifts. These massive galaxies could heavily influence the shape of the MS we observe, so it is vital to understand what is driving their infrared emission.

There could be several reasons why galaxies with quiescent colors have significant emission in the infrared, including (1) improper classification of star-forming galaxies possibly due to extreme dust obscuration, (2) elevated infrared luminosity from an AGN, (3) inaccurate absolute magnitudes due to catastrophic failures in photo-z's or low signal-to-noise photometry, or (4) "post-starburst" infrared glow due to dust heating from young stars (that is not related to the instantaneous star formation).

AGNs typically heat dust to very hot temperatures, so we expect any AGN contribution to infrared radiation to be predominantly in near- and mid-IR wavelengths, while far-infrared emission is likely due to star formation alone. Only about 10% of the IR-Q galaxies have been detected by Herschel, and the rest are 24 μm only detections, where AGNs may heavily influence the emission. However, only ∼1%–5% of the IR-Q galaxies are detected in the X-ray by Chandra, and only ∼2%–8% of the IR-Q galaxies have IRAC power-law colors indicative of AGNs, with significant overlap in those two populations, and the percentages are even lower when looking only at the 24 μm only sources. This suggests that radiation from an AGN is not fueling the infrared emission. The average SFR of the IR-Q galaxies with AGNs is 0.2–0.5 dex higher than the average SFR of all the IR-Q galaxies, so the presence of AGNs in these galaxies is likely just a reflection of the well-studied trend that AGN fraction increases with SFR or L_IR (e.g., Kartaltepe et al. 2010).

The rather high SFRs of the IR-Q galaxies suggest that they are indeed driven by star formation, and have been misclassified as quiescent. Man et al. (in preparation) stack the infrared emission from quiescent galaxies and find upper limits of SFR_IR  <  0.1–1 M_☉ yr⁻¹. The SFRs of the IR-Q galaxies in our sample tend to lie below the MS, but are all at least a factor of two to three higher than the upper limit from Man et al. (2014), which suggests that they are indeed still actively star forming. In addition, the infrared emission from IR-Qs is brighter than expected from a "post-starburst" glow (Hayward et al. 2014). It is unlikely that catastrophic photo-z errors or low signal-to-noise photometry are causing these misidentifications, as the sources have excellent photometry and well-constrained photometric redshifts (only 0.1% of the IR-Q galaxies have σ_{Δz/(1 + z)} > 0.15). Thus, the likely explanation for these sources is that they are actively star-forming galaxies that have been misclassified as "quiescent," and we include them in our analysis of the MS. We note, however, that the shape of the MS does not change significantly based on the inclusion or exclusion of these sources.

The relationship between SFR and M_* varies with stellar mass, with two distinct regions below and above the characteristic turnover mass, $\mathcal {M}_{0}$. What causes this change, and why does the turnover occur at about M_* ≈ 10¹⁰ M_☉ at all redshifts?

5.1.1. The Slope of the Main Sequence

5.1.2. Quenching in High Mass Galaxies

From our fits, we find strong evolution in $\mathcal {S}_{0}$, which parameterizes the overall scaling of the SFR/M_* MS with redshift.¹⁵ The scaling of the MS has been found in the literature to evolve as (1 + z)ⁿ, with the exponent n varying from 2.2 < n < 5 (Erb et al. 2006; Daddi et al. 2007; Damen et al. 2009; Dunne et al. 2009; Pannella et al. 2009; Karim et al. 2011). The value of n = 4.12 ± 0.10 we measure is among the steeper slopes seen in the literature.

It's thought that the redshift evolution of the main-sequence normalization is due, at least in part, to increasing gas content in galaxies at earlier cosmic times. However, measuring the gas content in galaxies can be difficult, especially in high redshift systems. Molecular hydrogen is notoriously difficult to detect, so many surveys instead probe the rotational transitions of CO and use locally calibrated CO-to-H₂ conversion factors, although this conversion factor may differ in starburst galaxies (Tacconi et al. 2008; Magdis et al. 2011; Magnelli et al. 2012). Another method to estimate gas content is to measure the dust mass from far-infrared or submillimeter photometry and convert to gas masses using an assumed gas-to-dust ratio (e.g., Magdis et al. 2012; Santini et al. 2014; Scoville et al. 2014). Magdis et al. (2012) find gas fraction evolves as (1 + z)^2.8, while recent ALMA observations suggest a steeper evolution of (1 + z)^5.9 (Scoville et al. 2014). Zahid et al. (2014) study the mass metallicity relationship and infer a much shallower evolution of gas mass M_g∝(1 + z)^1.35. Future studies will be needed to determine if the evolving normalization of the MS can be explained simply by an increasing gas supply in galaxies, or if other explanations such as increased merger rates or increased star formation efficiency are necessary to fully explain the observed evolution.