Star Formation at $4<z<6$ From the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH)

Using the first 50% of data collected for the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH) observations on the 1.8 deg$^2$ Cosmological Evolution Survey (COSMOS) we estimate the masses and star formation rates of 3398 $M_*>10^{10}M_\odot $ star-forming galaxies at $4<z<6$ with a substantial population up to $M_* \gtrsim 10^{11.5} M_\odot$. We find that the strong correlation between stellar mass and star formation rate seen at lower redshift (the"main sequence"of star-forming galaxies) extends to $z\sim6$. The observed relation and scatter is consistent with a continued increase in star formation rate at fixed mass in line with extrapolations from lower-redshift observations. It is difficult to explain this continued correlation, especially for the most massive systems, unless the most massive galaxies are forming stars near their Eddington-limited rate from their first collapse. Furthermore, we find no evidence for moderate quenching at higher masses, indicating quenching either has not occurred prior to $z \sim 6$ or else occurs rapidly, so that few galaxies are visible in transition between star-forming and quenched.


INTRODUCTION
Studies of star-forming galaxies at a wide range of cosmic epochs have revealed a strong correlation at fixed redshift between star formation rate (SFR) and stellar with α and β likely time-dependent. This relationship has been shown to hold over 4 -5 orders of magnitude in mass (Santini et al. 2009) and from z = 0 to z ∼ 6 (Noeske et al. 2007;Daddi et al. 2007;Karim et al. 2011;Kashino et al. 2013;Speagle et al. 2014). It is a tight correlation, with only ∼ .25 -.35 dex of scatter at any redshift (Daddi et al. 2007;Whitaker et al. 2012).
This paper extends the star-forming main sequence (SFMS) to massive galaxies at higher redshifts than previous studies by using the first data from the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH). This survey is obtaining 2475h (>6h/pointing) of Spitzer IRAC 3.6µm and 4.5µm observations on the two Hyper-Suprime-Cam (Takada 2010) ultra-deep fields COSMOS (Scoville et al. 2007;Sanders et al. 2007) and SXDS (Ueda et al. 2008). Here we use the first 50% of the data on the COSMOS field (∼ 4h/pointing) along with previously published 0.15-2.5µm data (Ilbert et al. 2013) to probe star forming galaxies with masses > 10 10 M ⊙ to z∼ 6. SPLASH is is a photometric survey that improves upon existing studies at 4 < z < 6. For these galaxies, the Lyman break allows quality photometric redshift determination, while SPLASH is complete and deep enough to greatly reduce biases that occur in studies specifically selecting for Lyman break galaxies (Lee et al. 2012;Speagle et al. 2014).
In § 2, we describe the SPLASH catalog, as well as the spectral energy distribution (SED) template fitting we use to produce photometric redshifts (photo-z's), stellar masses, and SFRs. § 3 describes the SPLASH view of the SFMS, including a lack of observed quenching at the high-mass end. Possible explanations for this lack of a turnoff are discussed in § 4, and in more detail in .

THE SPLASH DATASET
This work uses a revised version of the Subaru i + band selected Cosmological Evolution Survey (COS-MOS) catalog from Capak et al. (2007) to provide 0.15− 2.5µm photometry, with 10 6 spurious sources removed ), updated with intermediate band data, photo-z's, and physical parameters as described in Ilbert et al. (2010). The catalog was further augmented with significantly (1 magnitude) deeper z + band data taken with an updated Subaru-Suprime-Cam, Ultra-Vista (McCracken et al. 2012) imaging in Y, J, H, and K bands, and importantly the SPLASH IRAC 3.6µm and 4.5µm data. A full description of SPLASH is provided in Capak et al. (in preparation).
In this paper we use the first 50% of the SPLASH IRAC data, which covers a 1.2 degree diameter circle centered on the COSMOS field to 5σ ∼ 25.3 mag AB at both wavelengths. To reduce AGN contamination, X-ray detected sources (Brusa et al. 2010;Civano et al. 2012) were removed from the sample. To overcome source confusion (blending) we extracted the photometry using the i + band catalog as input for the IRACLEAN procedure described in Hsieh et al. (2012). The IRACLEAN photometry was compared with a T-FIT (Laidler et al. 2007) catalog in the CANDELS (Grogin et al. 2011;Koekemoer et al. 2011) area using an HST-WFC3 H band image as a prior and an EMPHOT (Conseil et al. 2011) catalog of the whole COSMOS field using the Ultra-Vista K band image as a prior. No systematic bias as a function of flux or position was found due to the photometric method.
The primary selection effect in our sample is the i-band catalog, which is used as a prior for the Spitzer photometry and as the basis for the photometric catalog. The use of this catalog limits us to z 6, when the Lyman break leaves the i-band. The depth limits us to unobscured star formation rates > 5 − 10M ⊙ /yr at z ∼ 4 − 6 (typical excitations are A V ∼ 0.4 mag for objects in our sample; cf. Bouwens et al. (2012a)). The typical depth of i-band selection is 26.5 mag (Capak et al. 2007;Ilbert et al. 2009). The second limit is the depth of the IRAC data, which introduces an age-dependent mass cutoff. For zero-age galaxies with a flat spectral energy distribution the IRAC data will limit the lowest measurable masses. For older galaxies, however, the depth of the Subaru i + image will limit the lowest measurable masses since they would not be in the i-band catalog. For stellar populations with approximately the age of the universe, our joint mass cutoff is noted in Figure 2. A third, more subtle effect is the method used to extract the IRAC photometry. If the photometry is not properly deblended, fainter galaxies will have systematically increased fluxes due to photometric crowding. The IRACLEAN algorithm corrects for this at least as well as other methods (TFIT, EM-PHOT), but we note that this may be affecting a small number of objects crowded by a source not in the i-band catalog.

Photometric Redshifts
We determine photo-z's using Le PHARE (Arnouts & Ilbert 2011) with the methodology described in Ilbert et al. (2013). Stellar masses were estimated with Bruzual & Charlot (2003) (BC03) models including strong emission lines with a mix of exponentiallydeclining (τ /100 Myr = [0.1, 0.3, 1, 3, 10, 30]) and delayed-τ (∆t peak /Gyr = [1, 3]) star formation histories (SFHs) with solar and half-solar metallicities at a range of ages spanning 0.05-13.5 Gyr to account for the effects of low metallicities and increasing SFHs on estimated physical parameters. Dust attenuation is modelled using the starburst Calzetti et al. (2000) curve and the λ −0.9 curve from Arnouts et al. (2013) Although the presence of a Lyman break should result in precise photo-z's, the 4000Å break and Lyman break can be confused resulting in catastrophic redshift failures. To check the quality of the photo-z's we crossreference them with spectroscopically confirmed sources (spec-z's), primarily from DEIMOS observations selected to be representative of z > 4 galaxies (Capak et al. 2010;Mallery et al. 2012). Out of 139 galaxies with specz > 4 and robust IRAC fluxes, 87 (63%) have successful photo-z determinations ( Fig. 1) to within 15%, with most of the rest mistakenly fit with at much lower redshift. Howvever, there are very few false positives in the SPLASH catalog. 6% of objects with photo-z > 4 and measured spec-zs lie at spec-z < 2 (indicating confusion between the two breaks), and after correcting for the significantly larger fraction of objects with measured spectroscopic redshifts at z < 2 (10% at z < 2 vs 4% at z> 4 with ch1< 25.3), this implies a contamination rate of 2%. So, the current photo-z analysis is preferentially scattering objects down from z ∼ 4 − 6 to z ∼ 0.5, (Fig. 1) with approximately 40% of SPLASH star-forming galaxies erroneously fit as low-redshift and removed from the sample. Tests using the two-point correlation and crosscorrelation functions confirm these fractions (Coupon et al., private communication).
We note that typical SFRs are far lower at z ∼ 0.5 than z > 4, so up-scattered objects fall well off the SFMS at z > 4 where they will not affect our further analysis.

Stellar Masses and Star Formation Rates
Once the best-fit photo-z's and extinctions are determined, Le PHARE fixes both and fits for a stellar mass and SFR based upon the BC03 templates with emission lines added at the best-fit redshift ). We compared these masses to those with several different SFHs, dust extinction laws, and emission line prescriptions similar to Arnouts et al. (2013) and found the results varied by less than 0.2 dex. Furthermore, we found consistent stellar mass estimates for objects at z > 4 when comparing these results with Ilbert et al. (2013), with the exception of heavily blended objects where our improved IRAC photometry had an effect. We note that although stellar mass estimates come from SED fits at all redshifts, at lower redshifts it is possible to measure SFRs independently. Although beyond the scope of this Letter, this is discussed extensively in Speagle et al. (2014). -A comparison of the best-fit photometric redshift and DEIMOS spectroscopic redshift for objects in SPLASH at z phot > 4 or zspec > 4. Due to confusion between the Balmer break and Lyman break, many objects with catastrophic errors are mistakenly excluded from the sample (scattered to z ∼ 0.5), but very few objects are scattered up, even though there are an order of magnitude more spectroscopic galaxies at low redshift.
At z ∼ 0 − 3, several studies have investigated SED-based SFRs and found there is generally good agreement between the SED and other SFR indicators (Salim et al. 2007;Arnouts et al. 2013;Carollo et al. 2013). To verify this at the higher-redshifts in SPLASH we tested full SED fitting-based SFR indicators against UV- (Meurer et al. 1999) and FIR- (Lee et al. 2013) based SFRs.
It was possible to derive reasonable UV SFRs using the IR excess -UV slope (IRX-β) relation (Meurer et al. 1999) for objects with detections in at least three bands in the rest-frame 1600-3000Å range. Galaxies with either negative A 1600 or unconstrained β were removed, leaving a sample of 404 galaxies. There is good agreement in the mean values of these UV and SEDbased SFRs over more than two orders of magnitude in SFR, but with a large scatter of 0.8 dex. This scatter is consistent with expectations at these high redshifts (Bouwens et al. 2013;Finkelstein et al. 2012) due to both the intrinsic properties of these galaxies and measurement error.
For galaxies with very high SFRs ( 500M ⊙ /yr) we use Herschel fluxes (Lee et al., in prep.) to estimate FIR-based SFRs and find generally good agreement. The scatter between the SED and FIR based SFRs is also large but consistent with other studies (e.g., Reddy et al. (2008)). We conclude that our SED based SFRs are reasonable, but may be under-estimating the actual scatter in the SFMS. One additional concern is that we only use one dust extinction law, but in reality there are likely many, and we cannot disambiguate at these redshifts.

Sample Selection
There are 2,002,871 sources in the SPLASH catalog, of which 1,319,197 are classified as stars or galaxies with well-constrained redshifts, masses, and SFRs. Of these, 7583 have a best-fit photo-z > 4 ( § 2.1). We removed any objects where quasar or stellar SEDs had better fits than the galaxy templates. Assuming conservatively that the stellar population has the same age as the universe, SPLASH is mass complete at log M * /M ⊙ > 10 at these redshifts, a limit that will be improved by over 1 mag with upcoming Hyper-Suprime-Cam (HSC) data in the visible.
To obtain a clean sample we chose a quality cut based on the SED fitting uncertainty, χ 2 /DOF < 5, rather than a single band signal-to-noise (S/N). In practice, this is primarily a selection for having many observed bands with S/N > 3. This leaves a sample of 3398 objects between redshifts 4 and 6 (2541/857 above/below z = 4.8, dividing the range into two bins of equal cosmic time). Limiting the sample to the mass-complete regime leaves 2152 objects (1513/639). We will use this smaller, mass-complete (sub)sample when deriving relationships between physical quantities to avoid biases. We include in our sample approximately 50 objects with no IRAC detection but otherwise well-constrained SED fits that pass the remainder of our cuts to avoid introducing a bias against particular SED shapes. Removing these objects has a negligible effect on the conclusions drawn in § 3.
Objects selected by the χ 2 cut are preferentially brighter in the optical, and may introduce a bias in the sample distribution. To test for this we did a Kolmogorov-Smirnov comparison between the inferred SFR and stellar mass of our χ 2 /DOF < 5 sample and the original sample. The two are consistent with being drawn from the same distribution.
As discussed in § 2.1, approximately 80 objects (∼ 2%) of our sample are expected to be low-redshift galaxies fit with z > 4 due to confusion between the 4000Å and Lyman breaks. At z > 4, it is expected that all galaxies are star-forming (Brammer et al. 2011;Behroozi et al. 2013), and z < 2 galaxies of any type scattered up to z > 4 will be assigned much lower SFR or even be selected as quiescent due to the shape of the typical galaxy SED. 73/3398 galaxies have SFRs below 10 M ⊙ /yr (70 with z < 4.8), while 46/3398 are found to be quiescent according to the Ilbert et al. (2013) criteria based on dust-corrected NUVrJ criteria. Because these two selections have little overlap and comprise a negligible fraction of the catalog, we choose to use all 3398 objects in our final sample.

THE STAR-FORMING MAIN SEQUENCE
The main result of this paper, the SFMS from SPLASH at 4 < z < 4.8 and 4.8 < z < 6 (bins of equal time interval), are shown in Fig. 2. In both redshift ranges, there is a strong correlation between stellar mass and SFR, qualitatively similar to that seen in approximately two dozen previous studies ) but with a different locus than at lower redshift. Almost all objects lie near this SFMS, with a small number of clear outliers at low SFR. Because of the limits of broad-band SED fitting, it is unclear without followup observations whether these outliers are candidates for quenched galaxies, catastrophic errors in photo-z's, and/or SED fitting, or have some alternative explanation.
The main sequence in each panel is divided by stellar mass into bins of width 0.1 dex, with the median for each bin shown in Fig. 3. For ease of comparison with other The "main sequence" for star-forming galaxies at photometric (a) 4 < z < 4.8 and (b) 4.8 < z < 6. A best-fit linear relationship is indicated by the blue dashed line in each panel. Mass and SFR completeness are indicated by the solid black lines. The magenta shaded region reflects an estimate of the increased sensitivity of SPLASH due to the addition of IRAC channels to existing multiwavelength data. Contours are drawn are equal intervals in number density. studies, medians above the mass completeness threshold that encompass > 25 objects are then fit with a power law (Eq. 1), with best-fit α and β indicated in Table 1. Using both fits, individual objects exhibit 1σ scatter of ∼ 0.24 dex about the main sequence, and agree well with Lyman-break selected objects from Lee et al. (2012) at similar/slightly lower redshifts and lower stellar masses.
Previous studies have disagreed on whether α is increasing (e.g., Santini et al. (2009)), constant (e.g., Karim et al. (2011)), or decreasing (e.g., Whitaker et al. (2012)) with redshift. We find a high-redshift slope of 0.78 ± 0.02 (Table 1), in good agreement with Lee et al. (2012) and the increasing-α prediction from the Speagle et al. (2014) literature compilation. The uncertainties given are statistical, but this good agreement may imply that the unknown systematic uncertainties are smaller than might otherwise have been expected.

High-Mass Galaxies and Quenching
As star formation is quenched, individual galaxies should turn off the SFMS, either lying at lower SFR or becoming quiescent and disappearing entirely due to selection. It has been previously suggested that highmass quenching is observed as a sub-linear SFR-M * relationship at high masses, but this may be due to selection effects (Noeske et al. 2007;Karim et al. 2011;Whitaker et al. 2012).
At the higher redshifts observed in SPLASH, the SFR main sequence is mass-and SFR-complete for the portion of the distribution indicated in Fig. 2. The SFR main sequence is reasonably approximated by a linear relationship for all masses above mass completeness and exhibits no evidence of a high-mass turnoff (Fig. 3). However, due to the depth of the current optical/NIR data, SPLASH is not sensitive to a population of fully quenched galaxies at this redshift. The star-forming galaxy main sequence median values (red dots) at (a) 4 < z < 4.8 and (b) 4.6 < z < 6, with contours from Fig. 2 superimposed. There is no decrease in star formation rate or any other evidence of quenching even for the highest-mass star forming galaxies. That is, the power-law form holds at all masses.

DISCUSSION
Several previous studies have examined the redshift evolution of the SFMS, finding that the slope α = d(log SF R)/d(log M * ) may decrease towards higher redshift. SPLASH allows us to check whether this trend continues to z ∼ 6. At lower redshift the evolution of α has been fit as a function of redshift (Whitaker et al. 2012) and of time . We find the time-based fit is better at high redshift as (1 + z) ∝ t 2/3 becomes more clearly non-linear (Table 1), which makes sense because the physical quantity affecting SFR is likely time, not redshift.
It is surprising that the SFMS at 4 < z < 6 has almost exactly the properties that one would extrapolate from lower-redshift behavior, with increasing specific SFRs (≡ ψ/M * ) for M > 10 11 M ⊙ galaxies. This trend cannot continue indefinitely because the virialization time for a galaxy depends upon its initial overdensity, and a typical galaxy with total baryonic mass of M > 10 11 M ⊙ (M ≥ M * ) will not virialize until z 10, with more massive galaxies having even lower virialization redshifts (Haiman & Loeb 1997). Furthermore, there is an effective Eddington limit on star formation for a given galaxy mass before star formation destroys the interstellar medium (Younger et al. 2008) of where the line-of-sight velocity dispersion in units of 400 km/s (σ 400 ), the interstellar medium opacity in units of 100 cm 2 /g (κ 100 ), and characteristic physical scale of the starburst in kpc (D kpc ) are likely of order unity. Thus, even under ideal conditions, it takes at least an additional ∼ 10 8 − 10 9 yr beyond virialization for the largest protogalaxies observed in SPLASH to attain M * > 10 11.5 M ⊙ . The universe at z ∼ 6 is slightly less than 10 9 yr old. At z ∼ 6 the formation of these galaxies is possible if they spend most of their time at near-maximum SFRs. Such galaxies would then not have a history of Best fit slope α and SFR at log M/M⊙ = 10.5 β for the star-forming main sequence in SPLASH. The expected fits according to the redshift and time dependence given by Whitaker et al. (2012) and Speagle et al. (2014), respectively, are included for comparison, as well as the best fit values for the Lyman-break selected samples of Lee et al. (2012) at z ∼ 3.9 and 5.0 corrected for extinction via Bouwens et al. (2012b).
sporadic starbursts, but rather near-continuous high-rate star formation from the moment of collapse until z = 6. This calculation, however, is only constrained to order of magnitude, so some variability in their SFHs is possible. However, even under ideal conditions this trend cannot continue far beyond z > 6, as there is insufficient time to build up stellar mass. Thus, we should expect that upcoming Hyper-Suprime-Cam (HSC) data, which is anticipated to extend the SPLASH catalog to z ∼ 8, should find a turnover redshift where the most massive galaxies are not observed, indicating we are seeing the rapid buildup of galaxies that have just completed their primordial formation and virialization. Otherwise, a substantially modified theory for early-universe structure formation will be required. HSC observations will also provide a substantial improvement in completeness in selecting high-mass, low-SFR galaxies. Such galaxies lying only slightly below the SFMS would be present in the current data, but improved Y-and z-band detection using HSC is necessary to select dusty and/or more quiescent high-mass, high-redshift galaxies. The current SPLASH catalog is sufficient to demonstrate that there is no population of turnoff galaxies lying just below the main sequence, but cannot determine whether this is because no galaxies have completed their star formation by z ∼ 6 or because these galaxies have rapidly changed their SED and are lost due to selection. Such galaxies have been observed to z ∼ 3.5 (Ilbert et al. 2013) and likely form in Eddington limited starbursts (Toft et al. 2014), but it is unclear when the quenching process began. Finding these first turnoff galaxies is essential to understanding when and how the most massive systems in the universe are quenched.
Taken together, our results provide strong hints that the most massive galaxies at high redshift all have a very similar history, one where they are forming stars at a nearly maximal rate from the first moment of virialization until a rapid turnoff at high redshift. This might turn out to be incompatible with the current consensus view that most galaxies form their stars through a series of individual starbursts, likely triggered by local events and environment and individual to the history of each galaxy. However, SPLASH is reporting only on the most massive galaxies, formed at the highest redshifts and very overdense initial fluctuations, so such galaxies are not necessarily typical. A natural next step is to follow up on this analysis, combining it with observational evidence about the formation process for galactic nuclei to try and build a consistent framework for galaxy evolution (cf. ), but such a model is beyond the scope of this Letter.