THE RELATION BETWEEN STAR FORMATION RATE AND STELLAR MASS FOR GALAXIES AT 3.5 ⩽ z ⩽ 6.5 IN CANDELS

This work uses multi-wavelength photometry from the CANDELS GOODS-S field (Grogin et al. 2011; Koekemoer et al. 2011). In addition to CANDELS, this work includes the Early Release Science (ERS), Hubble Ultra Deep Field (HUDF), and deep IRAC imaging in all four IRAC channels (3.6–8.0 μm) from the Spitzer Extended Deep Survey (Ashby et al. 2013) programs. Throughout, we denote magnitudes measured by Hubble Space Telescope (HST) passbands with the ACS F435W, F606W, F775W, F814W and F850LP as B₄₃₅, V₆₀₆, i₇₇₅, I₈₁₄, and z₈₅₀, and with the WFC3 F098M, F105W, F125W, F140W, and F160W, as Y₀₉₈, Y₁₀₅, J₁₂₅, JH₁₄₀, and H₁₆₀, respectively. Similarly, bandpasses acquired from ground-based observations include the CTIO/MOSAIC U-band; VLT/VIMOS U-band; the VLT/ISAAC K_s; and VLT/HAWK-I K_s.

We use fluxes from the catalog constructed by Guo et al. (2013). Guo et al. selected objects via SExtractor in dual-image mode with H-band as the detection image. As described in Guo et al., two versions of the catalog were constructed using SExtrator parameters that were (1) optimized in detection threshold and object deblending to identify faint, small galaxies (the "hot" catalog) and (2) optimized to keep large, resolved galaxies from being subdivided into multiple objects (the "cold" catalog). Both catalogs are then merged whereby any object in the "hot" catalog that falls within the isophote of a galaxy in the "cold" catalog is removed in favor of the "cold"-catalog object.

The HST bands were point spread function (PSF)-matched and the photometry is measured on the HST bands using the SExtractor double-image mode described above. For the ground-based and IRAC bands, the catalog uses TFIT (Laidler et al. 2007) to measure photometry of these lower-resolution images using the HST WFC3 imaging as a high-resolution template for the galaxies. We use the final version of the GOODS-S TFIT catalog which includes the new I₈₁₄ (CANDELS) and JH₁₄₀ (HUDF12) photometry.

In addition to the flux densities and uncertainties provided in this catalog, we include an additional uncertainty, defined to be 10% of the flux density of each object in band. This additional uncertainty accounts for any systematic uncertainty that may be related to the source fluxes themselves. This includes, for example, flat-field variations, PSF and aperture mismatching, and local background subtraction, many of which will be (to first order) proportional to the flux itself. The value of 10% was chosen such that the distribution of reduced χ² is ⩾1 and is justified based on arguments in Papovich et al. (2001). We add this uncertainty in quadrature to the measured uncertainties to estimate a total uncertainty on the flux density in each band for each object.

We use results from the recent CANDELS GOODS-S photometric-redshift project (Dahlen et al. 2013) which we briefly summarize here. A team of eleven investigators tested their individual photometric redshift fitting codes on blind control samples provided by the CANDELS team. A hierarchical Bayesian approach was then performed to combine the seven investigators' individual P(z) distributions to a final P(z) distribution for each object. The photometric-redshift (z_phot) is thereafter derived as the weighted mean of this distribution. Another sample was constructed as the median z_phot of all eleven individual results. The z_phot distributions from the medians and the combined P(z) methods both retained a lower scatter and outlier fraction than the results of any single investigator. Tests by Dahlen et al. (2013) showed that the hierarchical Bayesian z_phot method produces the best (smallest) scatter between the z_phot and spectroscopic redshifts. Finally, these methods were applied to the same CANDELS TFIT catalog (Guo et al. 2013) from which our data were obtained.

Figure 1 compares redshifts from the combined P(z) method with their highest-quality spectroscopic counterparts. The top panel exhibits a histogram of the number of objects used in our samples as a function of their photometric redshift. The bottom panel shows the ability of the photometric redshifts to recover known spectroscopic redshifts in the redshift range of this work.

Zoom In Zoom Out Reset image size

Unless otherwise specified, we use the median absolute deviation (MAD) to compute the equivalent standard deviation, σ_MAD, as the measure of scatter in given quantities (Beers et al. 1990), including the quoted scatter for redshift, stellar mass, and SFR. The σ_MAD is an analog for the 68% confidence, σ, if the error distribution were Gaussian and is therefore less sensitive to outliers (see Brammer et al. 2008). The MAD standard deviation in the photometric redshift accuracy ranges from σ_MAD/(1 + z) = 0.016 at z  ∼  4 to 0.028 at z  ∼  6, indicating that these photometric redshifts reliably recover known spectroscopic redshifts at high redshift.

Even in the highest quality spectroscopic redshift sample, there is a non-zero chance that some objects will have a misidentified z_spec due to a misinterpreted emission line or Lyman break (see discussion in Dahlen et al. 2013). So it is likely that some outliers are actually due to a misidentified z_spec rather than a poorly fit z_phot fit. The number of outliers where |z_spec − z_phot|/(1 + z_spec) > 0.1 are 2, 1, and 1 for z ∼ 4, 5, and 6, respectively (only 5%, 5%, and 11% of each sample). For the remainder of this work we use the z_phot catalog derived from the combined P(z) method, and substitute for high-quality (z_qual = 1) spectroscopic redshifts when available.

We selected objects according to their photometric-redshift (3.5 ⩽ z_phot ⩽ 6.5). This redshift range was chosen to be close to the redshift range of the traditional B, V, and i -drop samples. The lower redshift bound was chosen to avoid higher photometric redshift uncertainties, which may be due to a weaker Lyman break signal at z < 3.7 (see Dahlen et al. 2013, their Figure 11). Our samples have been cleaned from a total of 46 objects from X-ray (Xue et al. 2011), IR (Donley et al. 2012), and radio (Padovani et al. 2011) detected AGNs, as flagged by the Dahlen et al. (2013) photo-z catalog.

Objects with a best-fit SED with χ² > 50 are omitted from all samples. This cut removes objects with particularly poor fits, which comprise less than 4% of all objects. We interpret these objects as poor detections that do not represent the data well, and note that the removal of these objects do not impact the results of this work. The final sample includes 1728 objects with 3.5 < z < 4.5, 553 objects with 4.5 < z < 5.5, and 266 objects with 5.5 < z < 6.5, as illustrated in Figure 1. We refer to these as the z ∼ 4, 5, and 6 samples, respectively.

This work takes advantage of recent mock catalogs with synthetic photometry for galaxies from semi-analytic models (SAMs), as well as a semi-empirical dark matter and hydrodynamic simulation. We collectively refer to these as "the models." The benefit of comparing our derived results against these model galaxies is that the models incorporate realistic star formation histories and galaxy physics. Here we use these models for two comparisons. First, in Section 4.4 we derive stellar population parameters (SFRs and stellar masses) from the synthetic photometry for the model galaxies and compare to their "true" values as a test of our SED fitting procedures. Second, in Section 5.2.2 we use our derived SFRs and stellar masses from the CANDELS samples to compare to the models and interpret the SFR–mass relation and its scatter.

2.4.1. SAMs of Somerville et al. and Lu et al.

2.4.2. Semi-empirical Matching of Observed Galaxies to Dark Matter Halos of Behroozi et al.

2.4.3. Hydrodynamic Simulation of Davé et al.

One of the aims of this paper is to constrain the star formation history of the average population of galaxies at high redshift, z  >  3.5. Previous works have shown that broadband SED fitting offers no statistical preference between constant, rising, or declining star formation histories, even with broadband coverage spanning to the IR (Reddy et al. 2012). Furthermore, the star formation histories assumed in the templates can have non-negligible effects on the inferred SFRs, stellar masses, and ages (Lee et al. 2010). We addressed the shape of the star formation history as constrained by individual galaxies by running three separate fits using templates that assumed constant, declining and rising star formation history. The rising and declining star formation history templates ("τ" models) included a suite of e-folding times. We ultimately found no obvious χ² preference on the shape of the star formation history for individual galaxies.

Qualitatively, there has been some evidence to reject high-redshift star formation histories that decline with time. Previous studies have shown that high-redshift declining star formation histories would under-predict the sSFR at lower redshifts (Stark et al. 2009; González et al. 2010; Maraston et al. 2010). In addition, the instantaneous SFRs derived when assuming a declining star formation history will be under-produced by a factor of 5–10 as compared to direct estimates based off of UV-to-mid-infrared emission (Reddy et al. 2012). Other evidence against declining star formation histories comes independently from the SFR evolution of UV-luminous galaxies selected at fixed number density (Papovich et al. 2011). Finally, Pacifici et al. (2012) introduced a state-of-the-art SED fitting procedure with realistic, hierarchical mass-assembly histories and showed that declining τ-model histories do not well represent galaxies even at z < 2.

Although galaxies at z  >  3 likely have star formation histories that increase monotonically with time, we found it was impractical to use such models as the derived results are less physical. Our full justification for fitting individual galaxies with a constant star formation history is provided in Appendix A. Briefly, the BC03 stellar populations currently only allow for star formation histories that rise exponentially with time using simple parameterizations. At late times, such histories increase their SFR much faster than supported by observations. In Appendix A, we show that our modeling of synthetic photometry for galaxies from semi-analytic models recovers the most accurate stellar masses and SFRs when we adopt constant star formation histories. We interpret this as due to the fact that the SFRs in the models are approximately constant over the past ∼100 Myr (see, e.g., Finlator et al. 2006), and not consistent with exponentially increasing SFRs. Therefore, we fix the fitting templates to have a constant star formation history in our analysis of the CANDELS data for the remainder of this work. In a future work, we will explore possible improvements in parameters using models with star formation histories that increase as a power law in time (Ψ ∼ t^γ).

This section presents our method of incorporating nebular emission. Nebular emission is important because many galaxies at high redshift are observed to have intense star formation and high equivalent widths (EW) from emission lines (Erb et al. 2010; van der Wel et al. 2011; Atek et al. 2011; Brammer et al. 2013). Such strong nebular emission is able to enhance broadband flux by up to a factor of ∼2–3 in IRAC 3.6 and 4.5 μm bands (Shim et al. 2011). Previous studies have shown that the flux excess from high EW emission lines causes a systematic decrement in stellar mass and SFR inferred from SED fitting (Schaerer & de Barros 2009, 2010; Ono et al. 2010; Finkelstein et al. 2011; de Barros et al. 2014; Reddy et al. 2012; Stark et al. 2013).

3.2.1. Nebular Lines

The strength of a given emission line is dependent on the properties of both the gas cloud and the incident ionizing source. These properties include metallicity, ionization parameter, electron density, and number of ionizing photons. Inoue (2011) explored these parameters and the resulting strength of nebular emission in the regime of high-redshift galaxies by utilizing CLOUDY 08.00 (Ferland et al. 1998), which we use in our incorporation method.

After modeling a wide parameter space of seven metallicities, five ionization parameters and five hydrogen densities, Inoue (2011) reports 119 sets of metallicity-dependent emission line strength relative to Hβ. These line ratios, ranging from 1216 Å to 1 μm, are in close agreement with empirical metal line ratios (Anders & Fritze-v. Alvensleben 2003; Maiolino et al. 2008). We use the Inoue (2011) line ratios and include Paschen and Bracket series lines from Osterbrock & Ferland (2006) and Storey & Hummer (1995). Following Inoue (2011), we relate the Hβ line luminosity to the incident number of Lyman-continuum photons as

$$\begin{equation} {\rm L}_{{\rm H}\beta }\ =\ 4.78 \times 10^{-13} \frac{1-f_{\mathrm{esc}}}{1+0.6 f_{{\rm esc}}} {\rm N}_{{\rm LyC}}\ [{\rm erg\, s}^{-1}], \end{equation} \tag{ 1 }$$

where f_esc is the fraction of ionizing Lyman continuum (LyC) photons escaping the galaxy into the IGM and N_LyC is the production rate of hydrogen-ionizing photons (see also Osterbrock & Ferland 2006; Ono et al. 2010). The number of ionizing Lyman-continuum photons, N_LyC, depends on the age of the stellar population, and we take N_LyC from each BC03 SPS model for each age. It follows that 1 − f_esc is the fraction of LyC photons that ionize gas within the galaxy, which then produce the emission lines. The additional factor in the denominator of Equation (1) comes from a ratio of recombination coefficients (see Inoue 2011). Here, we equate the metallicity of the nebular gas to the metallicity of the SPS template (set as Z = 20% Z_☉ for all models, see above). Following the results of Erb et al. (2006b), we attenuate both nebular and stellar emission in the same manner (see Section 3.3 for details on attenuation).

The escape fraction has been measured to be low, i.e., f_esc ≈ 0 at low redshift z  ∼  1 (see Malkan et al. 2003; Siana et al. 2007, 2010; Bridge et al. 2010). At z ≳ 4, the IGM imparts a large optical depth to ionizing photons, making it difficult to constrain f_esc. Nestor et al. (2011) used z ∼ 3 Lyman break galaxies to study the high-redshift escape fraction, finding it to be consistent with f_esc ≈ 0.1. Finkelstein et al. (2012a) concluded that if galaxies are the main contributors to reionization, then the escape fraction must be f_esc < 0.34, or f_esc < 0.13 (2σ) at z ∼ 6 if the luminosity function extends to fainter galaxies than those observed, in order for the inferred ionization from galaxies to be consistent with the ionization background inferred from quasar spectra (Bolton & Haehnelt 2007). In addition, Jones et al. (2013) reinforced this claim by finding the covering fraction of neutral hydrogen in z ∼ 4 galaxies to be lower by 25% compared to z ∼ 2–3. From these results, it seems reasonable to assume a low, but non-zero, escape fraction at high redshift.

Here we consider two limiting cases. The first has f_esc = 1, for which all LyC photons escape the galaxy, preventing the creation of nebular emission and reverting the spectrum to the output of the SPS model. The second case has f_esc = 0, for which all LyC photons are absorbed and their energy is converted into the nebular emission spectrum. These two cases span the range of possibilities and allow us to study the effects of nebular emission on the inferred physical parameters. Nevertheless, given current constraints of f_esc ∼ 0.1 (see above), we expect the f_esc = 0 case will provide a more physical model for galaxies in our sample.

To illustrate the effect of nebular emission, Figure 2 displays four examples that include the best-fit SED models with and without nebular emission lines for galaxies. Depending on the redshift, emission from Hα and/or [O iii] will enhance the IRAC flux and can lead to highly different model-parameter values. The effect of nebular emission lines on the inferred stellar mass has a simple, qualitative explanation: the flux excess to the optical bands from nebular emission mimics a strong Balmer break that is typical for massive, older stellar populations. In this sense, when it is assumed that all of the observed broadband flux is produced by stars when much of it is produced by nebular emission, the inferred stellar mass will be overestimated.

Zoom In Zoom Out Reset image size

Similarly, Figure 3 illustrates how the specific emission lines ([O iii], Hβ, and Hα) affect the bandpass-averaged flux densities for the observed K_s and IRAC [3.6] and [4.5] bands from the best-fit SED models. The inclusion of [O iii], Hα, and Hβ lines raise the flux of these bands by up to a factor of approximately two. In Appendix B, we explore how the effects of nebular emission lines change the best-fit stellar masses and SFRs. However, as we show below, these changes are largely mitigated using our Bayesian formalism.

Zoom In Zoom Out Reset image size

3.2.2. Nebular Continuum

Recent work has suggested that the typically assumed Calzetti et al. (2000) attenuation prescription for star-forming galaxies is not ubiquitous (Reddy et al. 2012; Oesch et al. 2013; Chevallard et al. 2013; Kriek & Conroy 2013). The slope of the attenuation curve or presence of the UV dust bump at 2175 Å may be dependent on the galaxy type, geometry, metallicity, or inclination. However, galaxies at z > 4 currently lack sufficient observations to quantify these effects, so some attenuation prescription must be assumed. This work aims to test the effects of changing the type of assumed attenuation in order to gauge its impact on our broadband SED fitting procedure. In this subsection, we describe the two different attenuation prescriptions used in this work: the Calzetti et al. (2000) attenuation prescription ("starburst"-like attenuation hereafter), and the Pei (1992) attenuation prescription derived for the SMC ("SMC"-like attenuation hereafter).

Figure 4 shows four example best-fit SEDs of objects that emphasize the difference in SFRs for best-fit models using the SMC and starburst dust prescriptions. The starburst attenuation has a much "grayer" wavelength dependence in the UV than the SMC-like attenuation. This means the SMC-like attenuation curve has a much stronger attenuation at rest-frame, far-UV wavelengths λ ≲ 1200 Å, and a weaker attenuation across near-UV-to-near-infrared wavelengths λ ≳ 1200 Å, as shown in the top two panels of Figure 4. As stated above, bands shortward of Lyα are omitted in our procedure, so we do not fit where the difference between attenuation prescriptions is strongest.

Zoom In Zoom Out Reset image size

We find no obvious preference in χ² between the best-fit models for an SMC-like or starburst-like attenuation, and thus cannot as yet promote the use of one prescription over the other from this data set. However, we argue that the SMC-like attenuation could be invoked as a physical prior to reconcile the unphysical, extremely young stellar population ages that result from assuming a starburst-like attenuation. This method is preferred over, for example, increasing the minimum allowed age in the models (e.g., from ⩾10 Myr to ⩾60 Myr), which will not remove the preference of the fit to choose the youngest available age. A similar line of reasoning is used by Tilvi et al. (2013) to argue for SMC-like attenuation over starburst-like attenuation. Nevertheless, as we will show in Section 4.3, in our Bayesian formalism these differences arising from changes in the dust prescription are mitigated, and the dust-attenuation prescription has negligible impact on the results here.

Given a set of data for an individual galaxy which is a function of flux densities, D(f_ν), we derive the likelihood

$$\begin{equation} P(D| \Theta ^{\prime }) \propto \exp\, (-\chi ^2/2) \end{equation} \tag{ 2 }$$

where χ² is measured between the data, D, and a model in the usual way for a given set of model stellar population parameters, Θ' = (Θ{t_age, τ, A_V}, M_⋆). Note that the likelihood in Equation (2) is constructed based on linear fluxes. We then find the posterior probability density for any parameter given an observed set of data, D, and probability density using Bayes' theorem (see also Moustakas et al. 2013),

$$\begin{equation} P(\Theta ^{\prime }| D) =\ P(D|\Theta ^{\prime }) \times p(\Theta ^{\prime }) / \eta \end{equation} \tag{ 3 }$$

where Θ' represents the fitted parameters Θ and M_⋆, and η is a constant such that P(Θ'|D) will normalize to unity when integrated over all parameters (see Kauffmann et al. 2003). p(Θ') represents the priors on the model parameters, and is described further in Appendix C. As described in Section 3 and Table 1, we have adopted a prior (quasi-logarithmic) on the age from 10 Myr to the age of the universe for a galaxy's redshift, and we have adopted a prior (linear) that the attenuation is non-negative up to a maximum value. Further details on these priors and their effect on the fitting can be found in Appendix C.

We then derive posterior probability densities on individual parameters such as t_age, A_UV, etc. For example, the posterior on the age can be written as

$$\begin{equation} P(t_{{\rm age}}| D)\ = \int _{A_{{\rm UV}},\tau } P(t_{{\rm age}},A_{{\rm UV}},\tau | D)\ dA_{{\rm UV}}\ d\tau, \end{equation} \tag{ 4 }$$

where the integration is a marginalization over "nuisance" parameters, dust-attenuation, A_UV, and possible star formation histories/e-folding timescales, τ.¹⁷

The stellar mass must be treated differently because it is effectively a scale factor in the fitting process. In order to derive its posterior probability density we must integrate over all parameters, P(M_⋆|D)∝∫P(D|Θ')*p(Θ')dΘ. The mean and variance of the stellar mass can be computed as the first and second moments of the posterior. Similarly, the median stellar mass is defined as the value of M_⋆ such that the integral over the posterior from negative infinity to M_⋆ is equal to 50%, while the 68% confidence range can be calculated by integrating the posterior from the 16th to 84th percentiles.

We computed the posterior probability densities for all galaxies in our sample, including posteriors for the stellar mass, age, and attenuation, using the methods described above. Figure 5 shows examples of the cumulative posteriors on age, attenuation (or color excess, E(B − V)), and stellar mass for two galaxies in our sample: a relatively un-extincted and a relatively extincted galaxy. These objects are typical of those in the sample, where the posterior for age is typically broad, while that for the stellar mass is relatively tighter.

Zoom In Zoom Out Reset image size

In our analysis below, we will consider the relation between stellar mass and UV magnitude, as well as stellar mass and SFR for our full galaxy sample. Here, we discuss the relation between stellar mass and these quantities for an individual object, as it is illustrative. Figure 6 shows a two-dimensional probability density function between stellar mass and UV absolute magnitude. Here we take the M₁₅₀₀ from the conditional posterior on age and attenuation (similar to way we derive the posterior for stellar mass given the model parameters and data). Figure 6 also shows the posterior for M₁₅₀₀ and stellar mass individually.

Zoom In Zoom Out Reset image size

There is a weak covariance between M₁₅₀₀ and M_⋆, which results from the degeneracy in dust-attenuation and age. A galaxy with a redder rest-frame UV continuum has near equal likelihood to a model with an older, less extincted stellar population as to a model with younger, higher extinction. The two models produce a joint posterior that is anticorrelated between M₁₅₀₀ and stellar mass. The figure also shows the "main sequence" of the M₁₅₀₀–M_⋆ relation as derived from the full sample (see Figure 10). The joint posterior is approximately orthogonal to the main sequence, which implies that the likelihood scatter in each galaxy's M₁₅₀₀–M_⋆ plane will lead directly to scatter in the main sequence of the sample. We return to this point in the discussion of the SFR–M_⋆ relation below.

We derive the SFR from the model parameters in the following manner. We first determine the rest-frame UV luminosity at 1500 Å. At these redshifts, a large sample of detections redward of the rest-frame optical are unavailable, which can cause age and attenuation inferred from SED fitting to be degenerate quantities. These degeneracies can bias the M₁₅₀₀ inferred from the model parameters, especially if limited to best-fit values (see the discussion above and Appendix B). For this reason, we choose the closest observed band to rest-frame 1500 Å as a more observationally motivated value of M₁₅₀₀ because the broad-band photometry well samples the rest-frame UV portion of the SED. Galaxies have relatively blue rest-frame UV colors, so any corrections between the band closest to 1500 Å and the interpolated magnitude at 1500 Å are small (in the "extreme" examples of Figure 4 the differences are <0.1 mag). Furthermore, our tests show that none of our conclusions depend strongly on the manner we use to obtain M₁₅₀₀.

The 1500 Å luminosity is corrected for dust-attenuation using the median from the posterior of the attenuation at 1500 Å, or A_UV. The dust-corrected UV luminosity is converted to the SFR_UV using the ratio κ(t, τ) = SFR_UV/L₁₅₀₀. This is similar to the conversion given by Kennicutt (1998), but we account for variations in SFR_UV/L₁₅₀₀ owing to the age (t) and star formation history (τ) of the stellar population (see, Reddy et al. 2012). For each object, we use the median stellar population age from the posterior to calculate SFR_UV/L₁₅₀₀. However, we note that because most of these median ages from posteriors are >100 Myr for the galaxies in our sample and we have assumed constant star formation histories, in most cases κ(t, τ) is very similar to that of Kennicutt (1998).

We summarize the derivation of our SFR_UV mathematically as follows:

$$\begin{equation} \mathrm{SFR}_{{\rm UV}} = f_{{\rm CB}}\cdot \frac{4\pi D_L^2}{1+z} \cdot 10^{0.4\ A_{{\rm UV}}}\cdot \kappa (t,\tau) \end{equation} \tag{ 5 }$$

where f_CB is the flux of the closest band to rest-frame 1500 Å, D_L is the luminosity distance, A_UV is the median, marginalized attenuation at 1500 Å, and κ is the modified Kennicutt (1998) conversion that depends on age (t_age) and star formation history (τ).

The differences between the SFRs derived using Equation (5) and other common methods are described in Appendix B. In summary, methods that derive the SFR using the best-fit or direct UV luminosity slope exhibit higher scatter when compared to the marginalized SFR method from this work. This scatter stems from degeneracies between the young, dusty and old, dust-free solutions of a given SED, and photometric uncertainties (which affect the accuracy of measuring the UV spectral slope). We find the results of our method more robust as it reproduces SFRs from SAMs (see Section 4.4), and our method is relatively unaffected by model variations such as extinction prescription and/or nebular emission lines.

Figure 7 shows the SFR–stellar mass joint posterior for one object from our sample. As with the M₁₅₀₀–M_⋆ example, the covariance is roughly orthogonal to the star formation main-sequence, but there is more scatter because of the range in dust-attenuation (and, to a lesser extent, the stellar population age). The errors on the measured (extrinsic, or attenuated) M₁₅₀₀ are relatively small as they stem from photometric errors only, whereas the SFR depends on the UV luminosity corrected by the UV extinction, A_UV.

Zoom In Zoom Out Reset image size

SED models with higher A_UV have higher stellar-mass–to–light ratios (and therefore higher stellar masses at fixed UV luminosity). This induces some correlation in the SFR–stellar mass plane for each object. However, the covariance is mostly orthogonal to the expected direction of the SFR–stellar mass correlation, which implies that it contributes mostly to the scatter of the SFR–stellar mass relation and less to the correlation itself. In our analysis we take this covariance into account using Monte Carlo simulations (see Section 5.2.1).

Here we discuss our SED model assumptions and their impact on derived quantities such as SFR and stellar mass using our Bayesian method. In Appendix B, we show that these model choices have a significantly stronger impact on best-fit results, while the results using medians derived from posteriors are relatively unaffected.

The panels of Figure 8 show that the SFRs and stellar masses derived from the posteriors for the galaxies in our sample are rather insensitive to the choice of dust-attenuation prescription or the presence/exclusion of nebular emission (where we compare the results with f_esc = 0 or 1). In general, varying the assumed dust-attenuation prescription has a negligible impact on the derived SFRs and stellar masses (differences are <0.1 dex over the redshift range of our sample). Similarly, including emission lines has minimal effect on the SFRs (≲ 0.1 dex).

Zoom In Zoom Out Reset image size

There is some evidence that the stellar masses are reduced slightly when the models include nebular emission lines. This is in the same direction but weaker in magnitude as seen in comparisons of the best-fit models (e.g., see Appendix B). However, the effect is only a slight decrease of <0.25 dex, and is typically less than the measurement errors on mass derived from the posteriors. Therefore, the inclusion of nebular emission does not strongly impact the SFRs or stellar masses. For this reason, we have neglected exploring nebular emission over the full range of f_esc values in our analysis, and instead report results assuming f_esc = 0 (all ionizing radiation is absorbed and produces emission lines).

The results from Figure 8 can be directly compared to the results derived from best fits, which show stronger differences in these quantities when switching between the above model assumptions (see Appendix B). In contrast to using best-fit values, the Bayesian method uses the likelihood of all the models, and so even if there is a "highest likelihood" solution of low age, there are many good solutions with larger ages, and when marginalized, the latter dominate the posterior.

We tested the ability of our SED fitting procedure to reproduce accurately the SFRs and stellar masses in mock galaxy catalogs. We used synthetic galaxy photometry in the CANDELS bandpasses derived from the Somerville et al. SAM discussed in Section 2.4.1. The advantage of using synthetic catalogs from a SAM is that the models include realistic (and complex) star formation histories, as well as more sophisticated treatments of extinction (e.g., Charlot & Fall 2000). Therefore, the mock catalog from the SAM acts as a realistic observation where many of the model parameters (star formation history, extinction prescription) are known, and are distinct from the simpler models used in SED fitting to fit to the galaxy photometry. Comparing to the SAMs therefore provides a powerful test of our method to recover physical stellar population parameters, even when the physical details are unknown, such as the case of our observed galaxies. A similar but more comprehensive study was conducted among many CANDELS team members which compared the estimated and expected parameters from galaxies measured using SED fits (see B. Mobasher et al., in preparation).

The mock catalog from the SAM was filtered to a sample of 6000 simulated galaxies, evenly distributed across the mass and redshift range of our CANDELS sample. We took the synthetic photometry from the models and randomly perturbed them according to a Gaussian error distribution with a similar σ as the CANDELS data at a given band and magnitude. This process accounts for any systematic errors in creating the mock catalogs, and we process these fluxes the same way we process the data, even when fluxes are perturbed to negative values. These fluxes were used as input to the same SED-fitting procedure we applied to the real CANDELS samples. The masses and SFRs in the SAM were scaled to a Salpeter IMF, to match our procedure. We also fit to templates that exclude the effects of nebular emission (using f_esc = 1) because the SAMs also exclude nebular effects. For the test here, we show only the case where we fix the star formation history to be constant. Our tests (discussed in Appendix A) show that we recover the most accurate SFRs and stellar masses using constant star formation histories. Appendix A discusses how including additional star formation histories affects the SFRs and stellar masses.

Figure 9 compares the "true" stellar masses and SFRs from the SAM with those derived from using either the "best-fit" values or our method of taking the median, marginalized value from the posterior. This figure shows that taking advantage of the whole posterior with marginalized values produces less scatter in the recovering the "true" SAM values. This is because the maximum likelihood can be more sensitive to template assumptions than the median posterior values, as shown in Figure 8 and discussed in Appendix B. We note that there is a weak systematic where the fits slightly overestimate objects with low SFRs and slightly underestimate objects with high SFRs. The effect is mild, ranging by ±0.25 dex. This systematic could be due to differences in the extinction prescription and assumed star formation history between models used in the fit and those in the SAM.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 9 shows that the offsets in stellar mass and SFRs conspire to reproduce the accurate SFR–mass relation as expected from the SAMs. The parameters derived from our fits better reproduce the zero point, scatter, and slope of the relation than when using best-fit values. We also find no appreciable difference in the ability to recover the SFR–stellar mass relation across redshift. Our ability to test the success of our procedure, however, is limited to the maximum level of stochasticity in the star formation histories of the SAMs and we have not tested our procedure to observe its sensitivity to very bursty star formation histories. Nevertheless, these tests give us confidence that even in the presence of realistic photometric errors, we are able to derive a SFR–stellar mass relation from high-redshift galaxies in the CANDELS data that reproduces the intrinsic relation within these ∼0.25 dex uncertainties in SFR or stellar mass.

The panels of Figure 10 show the relation between the observed magnitude of the band closest to 1500 Å at the redshift of the galaxy (the rest-frame UV magnitude) and stellar mass at each redshift as two-dimensional histograms for our CANDELS sample. The median and σ of stellar mass in each bin of M₁₅₀₀ are given in Table 2. Here we show the observed UV absolute magnitude with no dust corrections in order to easily compare against previous studies (Stark et al. 2009, 2013; Lee et al. 2011, 2012). As discussed above, quantities used in this figure are median stellar masses derived from the marginalized PDF of each object (see Section 4 for details).

Zoom In Zoom Out Reset image size

We find a correlation between UV absolute magnitude and stellar mass, though this relation retains significant scatter. Recent evidence has suggested a relation with significant scatter between M₁₅₀₀ and stellar mass at high redshift, z ≳ 2 (Reddy et al. 2006, 2012; Daddi et al. 2007; Stark et al. 2009, 2013; González et al. 2011; Schaerer et al. 2013). The M₁₅₀₀–stellar mass trend in this work is weaker than the literature because we use an H₁₆₀-band selected catalog, which is closer to stellar mass than optically selected samples, as were used in the previous works listed above. This means that at fixed UV luminosity (or SFR) we are less sensitive to blue sources, which have higher mass-to-light ratios. Therefore, below our limiting stellar mass (10⁹/M_☉) we may be missing the bluer sources, as seen in Figure 10. At bright magnitudes (SFRs) our results agree with previous studies (that usually used z₈₅₀-band selected catalogs). It is at fainter magnitudes where our sources have fewer low-mass objects compared to the literature. Our results are also consistent with an independent analysis of the CANDELS catalogs, which used the same H₁₆₀-band selection (Duncan et al. 2014).

Regardless of the median relation, we consider the large scatter in M₁₅₀₀ at fixed stellar mass to mean one or both of the following. First, it could mean gas accretion is low such that galaxies undergo recurrent and stochastic star formation that leads to a range of M₁₅₀₀ at a fixed stellar mass (Lee et al. 2006, 2011). Second, galaxies at a fixed redshift and fixed stellar mass could exhibit a range of A_UV attenuations. In the second scenario, the observed scatter in the plane of M₁₅₀₀–M_⋆ would be largely diminished once we apply corrections for dust-attenuation. As discussed in the next subsection, the simulations favor the latter scenario.

The bottom right panel of Figure 10 illustrates an evolutionary cartoon depicting a model z ∼ 6 object of a given mass that is evolved forward ≈600 Myr to z ∼ 4. This is done by assuming an initial stellar mass (10⁸ M_☉), age (∼500 Myr, the average marginalized age of our entire sample), and zero dust (A_UV = 0). The strength and direction of three different star formation histories with the Bruzual & Charlot (2003) SPS code are shown: a constant SFR, a declining SFR (where Ψ ∼ exp (− t/τ), with τ = 1 Gyr), or a rising SFR (where Ψ ∼ t^γ using our results derived below from Section 6.2). As illustrated in Figure 10, only a rising star formation history naturally evolves galaxies along the median relation between stellar mass and M₁₅₀₀. Though this simple explanation does well to explain the UV-faint to UV-bright evolution, it offers little insight into the fate of the UV-bright galaxies at later epochs. It remains to be seen if some population of massive, UV-bright galaxies at z ∼ 6 quench their SFR such that we are missing a population of massive UV-faint galaxies at z ∼ 4.

Figure 11 shows the relation between the (dust-corrected) SFR and the stellar mass, where both parameters are derived from the fully marginalized probability density functions. Table 3 shows the median and σ scatter of log  SFR in bins of stellar mass from Figure 11. We measure a tight SFR–stellar mass relation (a "main sequence") for galaxies with log M_⋆/M_☉ >9, the mass completeness limit (e.g., Duncan et al. 2014). We explore how our SED fitting process could contribute to the correlation between SFR and stellar mass in Appendix C. This main-sequence in the SFR–mass relation has received much attention in the literature, and its existence implies that stellar mass and star formation both scale with the star formation history (Stark et al. 2009; González et al. 2010; Papovich et al. 2011). If true, then it follows that the gas accretion onto dark-matter halos at higher redshift is smooth when averaged over large timescales and stellar mass growth at high redshift is not driven by mergers (Cattaneo et al. 2011; Finlator et al. 2011). Our results support this picture.

Zoom In Zoom Out Reset image size

We fit a linear relation to the SFR–stellar mass relation as

$$\begin{equation} \log (\mathrm{SFR}/{M}_\odot \ \mathrm{yr}^{-1}) = a\ \log (M_\star /{M}_\odot) + b \end{equation} \tag{ 6 }$$

where a is the slope of the relation and b is a zero point. The fitted values for a and b are given in Table 4. We also show the fitted values for b when the slope is fixed to be a = 1, since the slope and intercept are often degenerate. We find that the slope and normalization in the SFR–mass relation shows no indication for evolution, with slopes of a =0.54  ±  0.16 at z  ∼  6 and 0.70  ±  0.21 at z ∼ 4. Furthermore, the scatter in SFR at fixed stellar mass shows no evidence for evolution, with a range of σ(log SFR/M_☉ yr⁻¹) = 0.2–0.4 dex from the median.

We must consider the possibility that the scatter in SFR at fixed mass is higher, and we are simply missing galaxies with low SFR due to incompleteness. We consider this unlikely because even if star formation ceased in some fraction of the galaxies, the galaxies would require 0.5–1 Gyr to have their SFR drop below a detectable threshold in the WFC3 IR data. These timescales are comparable to the period of time spanned by our subsamples (i.e., the lookback time between z = 4.5 and 3.5 is only 480 Myr), so it seems unlikely galaxies would "instantly" move from the observed SFR–mass sequence to undetectable values. For example, if such low-SFR objects existed at z = 4 their progenitors should be seen at z = 5 and 6 as they are fading, inducing a larger scatter in SFR–stellar mass. This work finds no evidence for such a population in our sample. Parenthetically, we note that some studies report evidence for massive, log M_⋆/M_☉ > 10.6 dex, quiescent galaxies at z ∼ 3–4, but this population lies at stellar masses above those in our sample (Straatman et al. 2014; Muzzin et al. 2013; Spitler et al. 2014).

We note that our SFR–stellar mass relation is tighter than our M₁₅₀₀–stellar mass relation (Figure 10), and we find that this can be explained by a correlation between stellar mass and our derived dust-attenuation (there is no correlation between derived attenuation and M₁₅₀₀). For example, objects at masses of 10^8.5, 10^9.5, and 10^10.5 M_☉ have median marginalized E(B − V) values of 0.05 ± 0.03, 0.13 ± 0.07, and 0.32 ± 0.18, respectively. This relation accounts for the differences in the scatter seen in Figures 10 and 11.

5.2.1. Constraints on the Intrinsic Scatter in the SFR–Mass Relation

Before comparing against models, it is necessary to understand how much of the scatter in the SFR–mass relation is intrinsic to the galaxy population and how much is a result of observational errors in SFR and stellar mass. To a simple approximation, the observed scatter (yellow in Figure 11) is a combination of the intrinsic (true) scatter and the measurement errors added in quadrature,

$$\begin{equation} \sigma _\mathrm{observed} =\left(\sigma _\mathrm{intrinsic}^2 + \sigma _\mathrm{errors}^2\right)^{1/2}. \end{equation} \tag{ 7 }$$

The SFR–mass joint probability density is broad, with covariance between the SFR and stellar mass (e.g., Figure 7). Because we calculate the posterior probability density functions for both the stellar mass and SFR for galaxies in our samples, we are able to estimate how correlations in these parameters contribute to the scatter and slope of the SFR and stellar mass relation. Here we use a Monte Carlo simulation to estimate σ_errors, and to determine how these errors affect the slope of the SFR–stellar mass relation.

We set up the Monte Carlo as follows. As discussed above, the SFR and stellar mass posteriors are covariant because both involve the dust-attenuation, A_UV, where models with higher A_UV have higher SFR from dust corrections, and higher mass-to-light ratios, which produce higher stellar masses. For each galaxy in each subsample at z = 4, 5 and 6, we randomly sample the galaxy's posterior density function of A_UV to find a new UV attenuation, A_{UV, i}. We then compute the conditional posterior for the stellar mass, P(M_⋆|A_{UV, i}). Next, we derive the SFR from Equation (5) and the medians of SFR in bins of stellar mass are re-calculated. This process is repeated 10⁴ times for each galaxy to generate 10⁴ new realizations of our galaxy sample. We then calculate at each stellar mass the median SFR and compute the σ_MAD from the distribution of medians. The scatter in log  SFR from this Monte Carlo is shown in Figure 11 and Table 3.

The σ_MAD scatter in the SFR from the Monte Carlo simulations is comparable to the observed SFR scatter, σ_observed, in most bins of mass and redshift. The scatter in the SFRs at fixed stellar mass from the Monte Carlo are shown in Figure 11 and given in Table 3. We make the approximation that σ_errors = σ_MAD from the Monte Carlo. We subtract these in quadrature from the observed SFR scatter to estimate the intrinsic scatter in SFR at fixed mass using Equation (7). We find that the average intrinsic scatter in SFR across the mass bins to be σ = 0.26 ± 0.04, 0.23 ± 0.10, and 0.34 ± 0.11 at z ∼ 4, 5, and 6, respectively. In some instances, the measurement errors from the Monte Carlo accounts for more scatter than the observed scatter, in which case there is no meaningful constraint on the intrinsic scatter. In this case, we take the Monte Carlo measurement scatter alone as a conservative limit on the intrinsic scatter (as some of the errors on the derived quantities must arise from the intrinsic scatter). This has implications for the gas accretion rate that we discuss below in Section 6.1.

The above test ignores the effects of our photometric redshift uncertainties since redshift is a fixed parameter during the fitting process. We constructed the following test to determine the effects of redshift uncertainties on SFR and stellar mass. We randomly selected 100 objects from each redshift sample and performed a Monte Carlo on their redshift uncertainty. In the Monte Carlo, each object's redshift was re-assigned according to a Gaussian error distribution with a σ equal to the object's 68% photometric redshift uncertainty. Then, we derived the stellar masses and SFRs in the same manner as with the data, fixing the redshifts to be the new redshift values. We calculated the medians in the SFR–stellar mass relation, as in as in Figure 11, for each of 10⁴ realizations of this process. Finally, we found the median and σ of SFR in each stellar mass bin from the distribution of SFR medians that each Monte Carlo realization produced.

Redshift errors produce a higher median log  SFR of <0.1 dex per stellar mass bin, and the redshift errors can contribute as much as ∼0.1 dex to the scatter in every stellar mass bin (usually it is much smaller, contributing <0.03 dex for 50% of the stellar mass bins). Therefore the redshift uncertainties do not significantly contribute to the error budget of the SFR–stellar mass relation.

We cannot rule out the possibility that a population of dusty, low-SFR galaxies are missing from our sample, which would attribute more scatter to the SFR–stellar mass relation. Indeed, some recent studies find evidence that such a population may exist at high redshift, at least at high stellar masses (Spitler et al. 2014; Man et al. 2014). Furthermore, at low stellar masses, our sample may be biased toward objects experiencing recent "bursts" of star formation (Schreiber et al. 2014). A deep investigation with ALMA is needed for further confirmation (Schaerer et al. 2014). However, our redshift range limits the data to rest-frame UV-to-optical wavelengths, and we defer the search for such a population for future work.

5.2.2. Comparison to the SFR–Stellar Mass Relation in Models

This subsection compares the results of the SFR–stellar mass relation in the previous section to results of recent SAM (Somerville et al., Lu et al.) semi-empirical dark matter abundance matching (Behroozi et al. 2013b), and hydrodynamic (Davé et al. 2013) simulations. Each of these simulations were briefly summarized in Section 2.4 and are collectively referred to as "the models."

Figure 12 shows the scatter in the SFR–stellar mass relation for each of the models as compared to the observed median and scatter in the SFR at fixed stellar mass (with data and errors bars identical to those in Figure 11). For each model, the median and scatter in the SFR at fixed mass was computed in the same way as the data (with all models converted to a Salpeter (1955) IMF, as assumed in this work).

Zoom In Zoom Out Reset image size

The SFR–mass relations from each model are in general agreement with each other and imply a tight relation between SFR and stellar mass exists for galaxies at high redshift. The SFR–mass relations from the models are also very similar to the observed relation derived from the data. Though some models predict a steeper slope of near unity (a  ≃  1), higher than measured in this work (a  ≃  0.6), the difference is negligible as it is within the errors. The observed offset in the zero point between the models and observations is likewise insignificant as the zero point is strongly anticorrelated with the slope in the linear fits of this work. In addition, from Figure 9 it was shown that the recovered SFRs from tests with the mock data tended to under- (over-) estimate the SFR of model galaxies with high (low) SFRs. No attempt was made to correct this for systematic. Therefore the similar offsets seen in Figure 12 between the high-SFR objects may imply that the data and SAMs are in even closer agreement.

As listed in Table 4, the scatter in the SFR at fixed mass is typically 0.1–0.2 dex in the models, whereas the limits on the intrinsic scatter from the data are σ(log SFR) < 0.2–0.3 dex (see Section 5.2.1). Therefore, both the observations and models support the conclusion that the SFR in galaxies at 3.5 < z < 6.5 at fixed mass (for log M_⋆/M_☉ > 9 dex) scales nearly linearly with increasing stellar mass and does not vary by more than a factor of order two. We explore the implication this has for the net gas accretion rate below in Section 6.1.

The fact that there is a tight SFR–stellar mass relation implies that the SFR scales almost linearly with stellar mass for the galaxies in our sample at 3.5 < z < 6.5. Because the SFR is (to a coarse approximation) the time derivative of the stellar mass, this implies that the SFR is an increasing function with time. We explore this further below in Section 6.2. Furthermore, the tightness of the scatter indicates that there is little variation in the SFR at fixed stellar mass. One caveat is that the SFRs are based on the galaxies' UV luminosity. Therefore, the SFRs that we measure are the "time-averaged" over the time it takes for the UV luminosity in galaxies to respond to changes in their instantaneous SFR. For the UV luminosity this timescale is approximately 30–100 Myr (e.g., Salim et al. 2009). Recent simulations have shown that the scatter is highly sensitive to the timescale of the SFR indicator (Hopkins et al. 2014; Domínguez et al. 2014). Therefore, the tightness in the SFR–mass relation of this work is conditional on the timescale associated with UV SFRs. With that in mind, the scatter observed in this work implies that galaxies at fixed mass in our sample have similar star formation histories when averaged over this timescale.

The scatter in the SFR–mass relation has important implications for net interplay between gas accretion into halos and galaxy feedback and outflows at these redshifts. The gas-accretion rate is a crucial piece of physics in galaxy formation models, and measures of the scatter in the SFR–mass relation are therefore an important test to constrain simulations and SAMs. As discussed above, the SFR is expected to track the net accretion rate/outflow rate of baryonic gas into the galaxies' dark-matter halos (Kereš et al. 2005; Dekel et al. 2009; Bouché et al. 2010; Lu et al. 2014). Because the SFR–mass relation is tight, it then follows that the gas accretion rate has little variation (σ  ∼  0.2–0.3 dex, or a factor of <2) at fixed stellar mass. Therefore, this favors a relatively smooth gas accretion process for galaxies at 3.5 < z < 6.5, at least above log M_⋆/M_☉ > 9 dex.

As discussed above, a tight, linear relation between SFR and stellar mass implies an SFR that increases with time. In this section, we study the SFR history directly to see if it is consistent with the observed SFR–stellar mass relation. This is achieved by tracking the evolution of the progenitors of z = 4 galaxies by selecting galaxies at different redshifts based on their number density.

Many studies have shown that a constant (comoving) number-density selection can trace the progenitor and descendant evolution both to relatively low and high redshifts (e.g., van Dokkum et al. 2010; Papovich et al. 2011; Lundgren et al. 2014; Leja et al. 2013b; Patel et al. 2013; Tal et al. 2014). In addition, recent studies have suggested using an evolving number-density selection to better track the progenitor populations of galaxies (e.g., Leja et al. 2013a; Behroozi et al. 2013a). Here, we use the parameterization of Behroozi et al. (2013a), who provide simple functions to track the number density evolution of the progenitors of galaxies.

This number density evolution is used to select the progenitors of galaxies in our sample. Figure 13 shows the cumulative stellar mass functions for the galaxies in our 3.5 < z < 6.5 CANDELS samples in bins of redshift. The results of Duncan et al. (2014) show that the objects in this field are complete for masses greater than log M_⋆/M_☉ = 8.55, 8.85, and 8.85 dex at z ∼ 4, 5, and 6 respectively. We assume a survey area of 170 arcmin² as described by Koekemoer et al. (2011) and the co-moving volume is calculated at each redshift assuming an uniform depth across each field. These cumulative functions are used to determine the stellar mass at which the galaxies of that redshift range achieve a given evolving number density as described by Behroozi et al. (2013a). As indicated in the figure, we take the galaxies with stellar mass log M_⋆/M_☉ = 10.2 dex at z = 3.75, which have a number density of log n/Mpc⁻³ = −4, and identify the galaxies at higher redshift that have a stellar mass that corresponds to the appropriate (de-evolved) number density at that redshift.

Zoom In Zoom Out Reset image size

Figure 14 illustrates our criteria to select objects according to an evolving number density. We use the Figure 13 stellar mass limits to find a best-fit curve across 3.5 <z < 6.5. Then, we select objects from our data that are ±0.25 dex in stellar mass about this relation. For the remainder of this work, we refer to these objects as "evolving number density–selected." We will later compare these briefly to objects selected at "constant number density" as such samples have received attention in the recent literature.

Zoom In Zoom Out Reset image size

Figure 15 shows the average SFR as a function of redshift for the evolving number density–selected galaxies. The SFR clearly increases as a function of time (decreasing redshift). We fit this evolution with a power law, Ψ(t) ∼ (t/τ)^γ where γ = 1.4 ± 0.1 and τ = 92 ± 14 Myr. If our sample is incomplete at low stellar masses (log M_⋆/M_☉ < 9.5 dex), then this would influence the measured power, γ. Based on Figure 13, the lower-mass objects will suffer greater incompleteness for objects of low SFR. This could mean that the intrinsic power-law slope is steeper than the one we measure here. We also note that we observe little difference between this evolution and that derived from using a constant-number density selection.

Zoom In Zoom Out Reset image size

Lastly, we can explore whether the SFR evolution derived above produces an average SFR–mass relation as measured from the data. We took the SFR history derived above and compute the resulting stellar mass and SFR for a stellar population starting at z = 6 and evolving to z = 4 using the BC03 SPS models (see discussion in Section 5.1 and bottom right panel of Figure 10). We plot the resulting SFR–mass relation in Figure 16 along with the medians derived above in Figure 11. Formally, the star formation history derived here corresponds to a galaxy with stellar mass log M_⋆/M_☉ = 10.2 dex at z = 3.75. Looking only at that point, the SFR–mass relation inferred from the SFR history matches the SFR–mass relation at z ∼ 4 remarkably well. (This is not circular because the stellar mass evolution is measured from the SFR evolution and not from the data itself.)

Zoom In Zoom Out Reset image size

Lastly, we use our derived SFRs and stellar masses to study the evolution of the sSFR. We first explore the evolution with redshift at fixed stellar mass as this has received attention in the literature (even though it is clear that galaxies with such high SFRs will not remain at the same stellar mass over this redshift range). We then consider the evolution in sSFR for the evolving number density-selected sample described above, as this will better track the evolution in galaxy progenitors across this redshift range.

Figure 17 shows the evolution of the sSFR for objects selected with constant stellar mass: within ±0.5 dex of log M_⋆/M_☉ = 9. The sSFR for objects at this stellar mass increases with increasing redshift, with high scatter. This evolution is consistent with hydrodynamic simulations (Davé et al. 2011b), SAMs (Somerville et al., Lu et al.), and other recent observational results from the literature (Stark et al. 2013; González et al. 2014).

Zoom In Zoom Out Reset image size

Figure 18 shows that the sSFR increases with redshift when galaxies are selected with an evolving number density. The evolution of the simple cosmological accretion models¹⁸ where sSFR ∼ (1 + z)^{ϕ = 2.5} (e.g., Neistein & Dekel 2008) is consistent with the sSFR evolution found in this work (Figure 18), which has ϕ = 3.4 ± 2.5 from z = 3.5 to z = 6.5. There is a weak difference between the evolution of the evolving number density-selected sample and the sample at constant mass in Figure 17; the sSFR evolution for the evolving number density–selection is shifted to lower values as a function of redshift. This is basically a confirmation of the fact that the SFR–mass relation produces an sSFR that is nearly independent of the stellar mass: sSFR = SFR/M ∼ M^{a − 1}, where sSFR ∼M^{−(≈ 0.4)} using our results in Table 4. This implies that samples selected at constant stellar mass or an evolving stellar mass (from a number density selection) return only slightly different sSFR values because the sSFR is a slowly changing function with stellar mass, as recently found by Kelson (2014). Consequently, this permits mass-independent modeling of the specific accretion rate, and therefore the sSFR (see Dekel & Mandelker 2014, for a detailed discussion).

Zoom In Zoom Out Reset image size

One way to explain the slight offset of the observed sSFR to lower values than the predictions from the SAMs is that some feedback mechanism may be present to hinder SFR from tracing halo or stellar mass growth. Gabor & Bournaud (2013) recently showed that cold streams of gas can increase the velocity dispersion in star-forming disks. They show that at z > 3 this increased turbulence causes the gas mass to stay higher, but reduces the SFR/gas mass by a factor of two. Assuming gas mass traces the baryonic-mass in galaxies at z > 3, then this could explain why the observed median sSFRs of this work are lower by about a factor of order two compared to other models.

The effects of including effects of nebular emission to the stellar templates are largely mitigated when using our Bayesian formalism (see Figure 8). In contrast, the nebular emission lines can strongly affect stellar masses and SFRs derived from best-fit model parameters. Figure 20 shows that in the redshift range of our galaxy sample, the presence (or absence) of an emission line in the ISAAC K_s, [3.6], and [4.5] bandpasses results in best-fit models that typically have lower SFRs and stellar masses by up to 0.5 dex (factor of three). This is because a galaxy with a redder rest-frame UV-optical color requires either an older stellar population or heavier dust-attenuation, with higher M_⋆/L ratios. Models where emission lines are present in the redder passbands reproduce the redder rest-frame UV-optical colors with lower stellar masses and SFRs. The effects of nebular emission are strongest in the redshift range that strong emission lines are present (see Figure 3). The median decrease is up 0.3 dex (factor of two) in both SFR and stellar mass for 5.5 < z < 6.5. As discussed in the literature, this can affect the interpretation of the galaxies (Schaerer & de Barros 2009, 2010; Schaerer et al. 2013; Ono et al. 2010; Finkelstein et al. 2011; Curtis-Lake et al. 2013; de Barros et al. 2014; Reddy et al. 2012; Stark et al. 2013; Tilvi et al. 2013).

Zoom In Zoom Out Reset image size

It is worth noting that strength of the nebular emission lines is highly uncertain for several reasons. One clear reason is that the model must use some escape fraction of ionizing photons, which is allowed to range between 0 ⩽ f_esc ⩽ 1. Another uncertainty is that in our parameterization, the strength of line emission is tied to the age of the model stellar population, and age is less constrained in SED fitting (see, e.g., Papovich et al. (2001) and Figure 5 above). Fits to galaxies using models with a starburst-like attenuation prescription (Calzetti et al. 2000) produce age distributions skewed heavily to the low ages (in some cases unphysically low as the ages are much less than a dynamical time; e.g., Yan et al. 2006; Eyles et al. 2007; Schaerer & de Barros 2009; de Barros et al. 2014; Reddy et al. 2012). Both the escape fraction and inferred ages are poorly known quantities which cause the effect of nebular emission lines to likewise be poorly constrained (see also Verma et al. 2007; Oesch et al. 2013). For these reasons it is advantageous to use SED fitting that is not strongly influenced by the presence or absence of such lines. As we show above (Section 4.3 and Figure 8), our Bayesian formalism is less affected by nebular emission lines in the models. Therefore, for the analysis in this paper we adopt models with f_esc = 0. In a future work, we will consider fully marginalizing over a range of escape fraction, although it seems unlikely to change the conclusions here.

Figure 21 shows the effects of the dust-attenuation prescription on the stellar masses and SFRs from the best-fit models. The choice of dust-attenuation prescription has a weak effect on stellar mass, where models using the SMC-like attenuation prescription have lower stellar masses by ∼0.1 dex in the median compared to the starburst-like dust-attenuation (although the spread is larger, up to ±0.2 dex).

Zoom In Zoom Out Reset image size

The choice of dust-attenuation prescription has a much stronger effect on the SFRs derived from best-fit models. Figure 21 shows that there is a strong trend in SFRs from the best-fit models: as the SFR increases, the models with starburst-like dust have significantly higher SFRs, by up to 0.5 dex, with a median of ∼0.3 dex (factor of two) compared to the best-fit solutions using an SMC-like dust. This is due to a combination of effects. First, there are high degeneracies between the inferred attenuation and age that arise from broadband SED fitting (discussed in the previous Appendix subsection). The assumed SFR depends on the stellar population age, especially at lower ages, where this leads to much higher ratios of in the SFR/L₁₅₀₀ ratio (see Appendix of Reddy et al. 2012) compared to the value typically assumed in the Kennicutt (1998) relation which assumes ages greater than 10⁸ yr.

Dust extinction and age are degenerate (negatively covariant) in the SED modeling and models with starburst dust typically have lower best-fit ages (Papovich et al. 2001). Therefore, the effects of higher dust-attenuation and higher SFR/L₁₅₀₀ ratio both contribute to a larger SFR for the case of starburst-like dust. The differences between starburst-like and SMC-like dust models are highest for models with highest SFRs, as shown in Figure 21. For objects with the highest SFR, the difference between the two prescriptions can be as high as ∼0.7 dex. Therefore, for best-fit models an assumed SMC-like attenuation causes starburst-attenuation-derived young, dusty, and high-SFR objects to be older, less dusty, and with lower SFR. This also reinforces the result of nebular emission in the Appendix above that simple template assumptions can significantly impact the best-fit SFR. Nevertheless, as we show above in Figure 8, in our Bayesian formalism these differences are mitigated, and the dust-attenuation prescription has negligible impact on the results here.

For these reasons, this work assumed an SMC-like attenuation for all objects in our sample. If we otherwise had very reliable age estimates we might assume, as Reddy et al. (2012), an "age"-dependent attenuation prescription, such that younger galaxies have SMC and older have starburst attenuation. However, the effects of nebular emission and assumed attenuation are strongest when adopting best-fit values. In contrast, we mitigate these effects by using the results from our Bayesian analysis, where the effects of changing dust-attenuation prescription are minimized (see Section 4).

The method used in the work to derive UV SFRs for high redshift galaxies is different from the typical methods found in the literature. The common methods include using a dust correction based on the UV spectral slope, f_λ ∼ λ^β, (Meurer et al. 1999; Madau et al. 1998; Kennicutt 1998) and using the instantaneous SFR from the best-fit stellar population model. As shown in Figure 22, both of these alternative methods show high scatter when compared to the marginalized SFR from the method of this work.

Zoom In Zoom Out Reset image size

The large scatter in Figure 22 is due to the fact that the scale of SFR is predominantly dependent on the treatment of the dust correction to UV luminosity, which is an unconstrained quantity in traditional SED fitting methods at high redshift. The Bayesian approach has the advantage in producing realistic ages, thereby reducing degeneracies with dust corrections. The median age for the SAM mock catalogs is log (t_age) = 8.48 ± 0.22 dex which resembles the distribution of marginalized ages found in this work for observed galaxies, log (t_age) = 8.73 ± 0.14 dex. Conversely, the distribution of best-fit ages is typically very dissimilar from the SAM and simulation predictions. This is because best fits typically find lowest χ² at the extreme end of parameter space (youngest ages, highest extinction). As a result, the median best-fit ages are lower with higher scatter log (t_age) = 8.06 ± 0.94 dex. This scatter results from the degeneracies between the young, dusty and old, dust-free solutions of a given SED and photometric uncertainties (which affect the accuracy of measuring the UV spectral slope, β), thus producing a wide range of acceptable SFRs.

In summary, when marginalizing over other nuisance parameters, the posterior on age returns more physical ages on the order of ∼350–750 Myr, reducing degeneracies in the derived dust corrections and thereby reducing the uncertainty in SFR. In addition, the method used in this work reproduces SFRs from semi-analytic models (see Section 4.4), and is relatively unaffected by model variations such as extinction prescription and/or nebular emission lines. This is additional evidence to favor the Bayesian approach to derive physical properties.