TRACING THE MASS GROWTH AND STAR FORMATION RATE EVOLUTION OF MASSIVE GALAXIES FROM Z ∼ 6 TO Z ∼ 1 IN THE HUBBLE ULTRA-DEEP FIELD

The galaxies we examine have been selected from new H₁₆₀-band imaging from the ultra-deep, near-IR WFC3/IR observations in the main HUDF field, part of the 192-orbit HUDF09 program (PI Illingworth: G0 11563). The HUDF09 program includes WFC3/IR images in three filter bands: Y₁₀₅, J₁₂₅, and H₁₆₀. The resulting images reach limiting magnitudes of ∼29 mag (5σ with 035 diameter apertures; see Bouwens et al. 2010a) with a native pixel scale of 006 pixel⁻¹ and a point-spread function (PSF) full width at half-maximum (FWHM) of ∼016. While deep imaging has also been obtained as part of the HUDF09 program for the flanking deep fields (HUDF09-1,2) within the larger HUDF05 (PI: Stiavelli), for the purpose of this analysis we incorporate only data from the primary HUDF09 field, which total 4.7 arcmin².

Imaging of the HUDF obtained with the ACS instrument on HST includes observations in four optical filter bands: F435W, F606W, F775W, and F850LP. The observations were obtained in a total of 400 orbits between 2003 September and 2004 January (Beckwith et al. 2006). The co-added images, now publicly available, have been matched to the pixel scale of the new WFC3/IR data. The respective limiting magnitudes of the B₄₃₅, V₆₀₆, i₇₇₅, and z₈₅₀ images are 29.7, 30.1, 29.9, and 29.4, as determined at 5σ within a 035 diameter aperture.

We include deep U-band imaging taken with the VIMOS instrument on the Very Large Telescope's (VLT's) Melipal Unit Telescope at Cerro Paranal Observatory in Chile. These data were obtained by the European Southern Observatory large program 168.A-0485 (PI C. Cesarsky), as part of the Great Observatories Origins Deep Survey (GOODS). The individual frames contributing to the co-added image were obtained in service mode between 2004 August and 2006 October. The combined U-band image reaches a 5σ limiting magnitude of 28.6, measured in a 12 diameter aperture with an average PSF FWHM of ∼08 (Nonino et al. 2009).

The K_s-band photometry in this work has been extracted from a combination of images taken with the Infrared Spectrometer and Array Camera (ISAAC) mounted on the Antu Unit Telescope of the VLT and the PANIC near-IR camera on the 6.5 m Magellan telescopes. The VLT/ISAAC exposures were taken in two programs (ESO program 73.A-0764; PI I. Labbe), one consisting of 8 hr on the GOODS Chandra Deep Field South and another 15 hr exposure in the UDF. The Magellan PANIC data include 15 hr of exposure. All K_s-band observations were reduced using standard procedures (Labbé et al. 2003). These observations have been combined and matched to the native pixel scale of the WFC3/IR images. The final combined image reaches a limiting magnitude of 27.4, with an average PSF FWHM of ∼036.

We include the super-deep observations from the GOODS-S Field, which were obtained using the Infrared Array Camera (IRAC; Fazio et al. 2004) on the Spitzer Space Telescope taken as part of the Spitzer Legacy Program. These data were taken with two pointings, each with an exposure time of approximately 23.3 hr. The ∼40 arcmin² of overlapping area with the HUDF received twice the exposure. We include imaging from channels 1 and 2, which are centered on wavelengths of 3.6 μm and 4.5 μm, respectively, each with a 1σ limiting depth of ∼27.

The PSF varies significantly between the images included in this analysis, such that the PSF FWHM generally increases with increasing central wavelength of the broadband imaging. If left uncorrected, this variance would produce a wavelength-dependent bias in our aperture photometry.

To alleviate this problem in the HST data, we degrade each image to match the PSF of the image with the broadest PSF, which in this case is the J₁₂₅ image. We begin by determining the PSF of the ACS and WFC3/IR images using a stack of seven non-saturated stars found in the field. Any detected neighboring objects are masked to reduce background fluctuations in the stacked image. We produce kernels for the stellar composite in each band using the lucy algorithm in IRAF and each image is convolved to match the J₁₂₅ PSF using the IRAF fconvolve routine. The FWHMs of the PSFs in the resulting images match to within 2%, as does the fraction of the total flux enclosed in a 10 pixel radius aperture.

Images with PSFs that are significantly broader than that of the J₁₂₅ detection image require more complicated PSF matching techniques. For the VIMOS, K_s band, and IRAC images, we therefore apply a different method, which allows for precise aperture photometry in cases where broad-winged PSFs can produce considerable contamination in the flux measurements of neighboring sources. We employ the same source-fitting algorithm described in detail in Labbé et al. (2006), Wuyts et al. (2008), and Whitaker et al. (2011). This method produces a model PSF for the image with the broader native PSF, which is then used to estimate the flux distribution of each source identified in the detection image segmentation map output by SExtractor, as described in the following section. For each individual object, the flux from neighboring sources is modeled and subtracted, allowing for a reliable aperture flux measurement of individual objects, even in reasonably crowded areas of the field. An example illustrating this PSF-matching technique, which we use to directly compare VIMOS, K_s band, and IRAC images with the H₁₆₀ detection image, is shown in Figure 1.

Zoom In Zoom Out Reset image size

In total, we measure K-band fluxes with a signal-to-noise ratio (S/N) >2 for 728 objects and 3.6 μm fluxes with an S/N >2 for 326 objects. For a small number of cases, particularly for objects in close proximity to very bright IRAC sources, blending problems are so severe that this procedure is unable to produce reasonable estimates of the flux in faint nearby objects. By visual inspection, 34 such cases were identified in the photometric catalog. Where present, this deblending failure generally produces an overestimate of the flux in the K_s band and IRAC bands. Since the space-based optical and near-IR photometry otherwise contains negligible contamination for most of these objects, they are retained in the catalog, but their fluxes in the K_s and IRAC bands have been set to –99. While these objects are retained in the catalog accompanying this work, we do not include them in the scientific analysis described in Section 5, as the K_s and IRAC flux measurements are critical to the accuracy of stellar mass estimates. The objects that we have removed from our analysis because of catastrophic deblending problems in the IR represent just 2% of the total data set and are drawn from a random distribution in redshift. Thus, we assume that any correction to our analysis due to their removal is negligible.

We use the SExtractor code (Bertin & Arnouts 1996) in dual-image mode to extract sources detected in the WFC3/IR H₁₆₀ image. In order to minimize false detections while maximizing the efficiency of faint source detection, we require a minimum detection area of 5 pixels with a threshold of 5σ significance above the background for both detection and analysis. We apply 64 deblending sub-thresholds with a minimum contrast parameter of 0.002. The background is estimated in GLOBAL mode, using a mesh size of 100 pixels and a filter size of 3 pixels.

The initial catalog output from SExtractor contains 1645 objects, which is reduced to 1553 after removing objects with coverage of <30% of the maximum in the H₁₆₀ weight map. This cut effectively removes spurious detections near the edge of the WFC3/IR imaging. An additional eight objects are removed after being identified as stars or artifacts.

In order to maximize the average S/N, circular apertures with a 10 pixel radius, corresponding to 06 at the pixel scale of the WFC3/IR images, are used to measure color fluxes for all objects. Total flux measurements are estimated in the H₁₆₀ detection image, using the SExtractor AUTO photometry and the following criteria:

$$\begin{equation} r_{{\rm aper}}= \left\lbrace \begin{array}{@{}rl}r_{{\rm Kron}}, &{r_{{\rm Kron}}>r_{{\rm color}}} \\ r_{{\rm color}}, &{r_{{\rm Kron}}\le r_{{\rm color}}} \end{array} \right. \end{equation} \tag{ 1 }$$

$$\begin{equation} H_{{\rm total}} = H_{{\rm AUTO}} \times \frac{1}{F_{r<r_{{\rm aper}}}}, \end{equation} \tag{ 2 }$$

where r_color = 06, r_Kron is the radius of the SExtractor circularized Kron aperture in units of arcseconds, and $F_{r<r_{{\rm aper}}}$ is the aperture correction used to recover flux falling outside the total aperture radius, r_aper. This correction is determined by the relative flux in the H₁₆₀ growth curve at the position of the circularized Kron radius of the AUTO aperture. Measurements of the total flux in each of the additional HST imaging bands are calculated by scaling the 10 pixel aperture fluxes in the following manner:

$$\begin{equation} F_{{\rm total}} = F_{{\rm aper}}\times \frac{H_{{\rm total}}}{H_{{\rm aper}}}. \end{equation} \tag{ 3 }$$

We derive errors for the space-based photometry by measuring the 1σ variance in a sampling of 2000 randomly placed circular apertures with a 10 pixel radius in each broadband image. Before placing the apertures, we subtract all detected objects using the SExtractor segmentation map as a mask. From the average error returned by this method, we determine the total error on the measurement of F_total for each object, scaling up the estimated 1σ color aperture error by the same factor of H_total/H_aper used to calculate F_total in Equation (3). Errors on the photometry in the U, K_s, and IRAC bands are determined by a similar sampling of the residual source-subtracted model images (e.g., the bottom-center panel in Figure 1).

For completeness, we include the catalog in the electronic version of the paper (see Appendix A).

As the HUDF lacks a sufficiently large spectroscopic sample representative of the range of galaxies contained in the field, neural network redshift codes are poorly suited to our analysis. Instead, we estimate photometric redshift using the EAZY code (Brammer et al. 2008) for all galaxies that have at least seven broadband flux measurements. EAZY operates by fitting a linear combination of template SEDs to each galaxy. The template set, which is described in detail in (Brammer et al. 2008), contains seven SEDs spanning a broad range in galaxy spectral types. The selection of these templates has been optimized to fit a wide range of galaxy properties, while minimizing degeneracies in color and redshift. In addition, an H₁₆₀-band flux prior is applied to limit catastrophic failures from remaining degeneracies.

For each object, EAZY produces a redshift probability distribution, which can be used to determine the photometric redshift in multiple ways. The discrete peak of this distribution is returned as z_peak and z_MC provides a redshift from Monte Carlo sampling of the probability distribution. For objects in which the redshift probability distribution is narrowly defined, the difference between these estimates is negligible. However, in cases where the probability distribution is broad, z_MC provides a more honest sampling within the error in the best-fit photometric redshift. For comparison, we show the photometric redshift distributions of z_peak and z_MC for the 1553 galaxies in the HUDF in Figure 2. While the results do slightly differ, a Kolmogorov–Smirnov test of the two distributions returns a statistic of 0.027, indicating that these differences may be considered to be insignificant for the average object in this sample. A comparison of the relative differences of the z_peak and z_MC estimates for objects matched to the FIREWORKS catalog implies that the difference between these estimates is negligible (differing on the 1% level). For the remainder of this work, we adopt z_peak as the photometric redshift estimate of choice, since the z_MC values are, by definition, impossible to exactly reproduce.

Zoom In Zoom Out Reset image size

In order to test the accuracy of the z_peak photometric redshifts, we cross match our catalog with the FIREWORKS catalog of Wuyts et al. (2008), the Grism ACS Program for Extragalactic Science (GRAPES; Pirzkal et al. 2004), and the Probing Evolution and Reionization Spectroscopically (PEARS; Straughn et al. 2008) catalogs (Rhoads et al. 2009). Among the 1553 galaxies in this H₁₆₀-selected HUDF catalog, 63 have spectroscopic redshifts in the aforementioned data sets, extending to z ∼ 5.5. The comparison of photometric and spectroscopic redshifts for the objects matching existing catalogs is shown in Figure 3. The mean difference calculated for (z_spec–z_peak)/(1+z_spec) over the range available is 0.017. We find a normalized median absolute deviation, σ_NMAD, of 0.055, as defined by

$$\begin{equation} \sigma _{{\rm NMAD}} = 1.48 \times \mathrm{median} \left \vert \left (\frac{\Delta z - \mathrm{median} (\Delta z)}{1+z_{{\rm spec}}}\right) \right \vert. \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

A total of 82 unblended objects in the H₁₆₀-selected galaxy catalog had photometric matches in the FIREWORKS catalog. A comparison of the photometric redshift estimates for these objects is also shown in Figure 3, with σ_NMAD = 0.048. The errors in the photometric redshifts we measure are comparable to measurements at z < 4 in other deep fields (σ_NMAD = 0.046; Brammer et al. 2008).

For galaxies with z > 3, the σ_NMAD with respect to the spectroscopic redshifts nearly doubles to σ_NMAD = 0.095 and the mean difference between spectroscopic and photometric redshifts increases to (z_spec − z_peak)/(1 + z_spec) = 0.07, indicating a potentially significant bias in the photometric estimates at higher redshifts. The scale of this worsening accuracy with redshift is consistent with the commissioning tests of Brammer et al. (2008), although the performance of EAZY beyond z = 4 has been limited to only a few galaxies. We have examined the best-fit SEDs for these high-redshift sources and found that the worsening agreement of the photometric redshifts with the spectroscopic measurements at z > 3 results from the difficulty of fitting the predominantly steep SEDs of young highly star-forming galaxies at these redshifts. While the trend worsens with redshift, the redshift bin sizes we choose for the analysis presented in Section 4 are sufficiently large (1 ⩽ δz ⩽ 2) to make these uncertainties negligible.

SPS models provide a means of estimating the evolution of an integrated galaxy spectrum, given an initial mass function (IMF), star formation history, and metallicity. This technique, first employed by Tinsley (1968), has shown dramatic improvement in recent years and is now a well-accepted method for deriving the estimated physical properties of galaxies from precisely measured SEDs.

The galaxy properties output by SPS modeling are inherently dependent on the input parameters, such as the IMF. While Chabrier (2003) and Salpeter (1955) IMFs produce similar colors for galaxies at z > 2, the Chabrier (2003) IMF predicts stellar masses that are ∼1.6 times smaller. The Chabrier (2003) IMF has been favored as of late and used in stellar mass estimates of high-redshift galaxies (e.g., Papovich et al. 2011). However, in a recent study of the stellar absorption features in local elliptical galaxies, van Dokkum & Conroy (2010) revealed an abundance of low-mass stars, contributing to more than 60% of the total stellar mass. This finding suggests that the IMF of massive star-forming galaxies in the early universe likely resembled a "bottom-heavy" IMF, more akin to, but steeper than, a Salpeter IMF.

We fit stellar population models to the measured SEDs at the redshift determined by EAZY using the FAST algorithm described in Kriek et al. (2009). We choose to employ a Salpeter (1955) IMF and the SPS models of Bruzual & Charlot (2003), as these seem to be most effective at modeling massive star-forming galaxies at high redshift. We also assume an exponentially declining (τ model) star formation history and solar metallicity. In estimating the SFRs, we apply the prescriptions of Wuyts et al. (2011), by requiring a minimum e-folding time of log(τ_min) =8.5, which has been shown to best reproduce the low-to-intermediate SFRs within the Bruzual & Charlot (2003) framework. For the estimation of stellar masses, we allow more freedom in the best-fit SPS models, setting log(τ_min) =7. We then estimate the specific star formation rate (sSFR) of each object using the combination of the results from these methods.

For each object, FAST determines the best fit to a six-dimensional cube of SEDs generated with a range in each of the following stellar population properties—age, star formation timescale, dust content, metallicity, and redshift—although we have fixed the redshift and metallicity in each case. In addition to the best-fitting results, FAST computes 68% confidence intervals on each stellar population parameter. These are determined using Monte Carlo simulations, which re-fit each observed SED 100 times after perturbing the flux measurements randomly within the photometric errors. As a result, we recover realistic random errors for each galaxy's estimated stellar mass, SFR, and age, with which we may now investigate their cosmologically averaged evolution in the following section.

We note that we have required a fixed metallicity, which should be a poor assumption, particularly for low-mass galaxies at high redshifts. Fitting the sample with metallicity as a free parameter indicates that a fixed solar metallicity will affect the best-fit estimated SFRs of the galaxies as a function of redshift. However, while this bias is clearly evident within the full photometric catalog, the galaxies above the lower mass limit of our analysis in the following section (log M_*[M_☉] > 9.5) show no significant departure from the results fit with a fixed solar metallicity. Thus, the influence of evolving metallicity should not inflict a substantial bias on the analysis presented in this work.

In order to best approximate the evolution of similar galaxies across a wide range in redshift, we have chosen to divide the galaxies into samples of constant cumulative number density, using a similar method to that utilized by Conroy & Wechsler (2009), van Dokkum et al. (2010), and Papovich et al. (2011). Contrary to galaxies selected by parameters that strongly evolve with time, such as mass, color, or luminosity, number density-selected samples more closely track the cosmologically averaged progenitors and descendants of galaxies with minimal bias over a range in redshift. This allows for the measurement of the average growth in stellar and gas mass and average SFRs for galaxies at a non-evolving number density over time.

As shown by van Dokkum et al. (2010) and Leja et al. (2013), even in the presence of mergers, samples selected by this method are relatively pure tracers of galaxy evolution. To precisely quantify the errors in cumulative number density selection caused by merger activity, Papovich et al. (2011) applied this selection to the simulated merger trees produced by the Millenium Simulation (Springel et al. 2005). Their results indicated that a constant cumulative number density selection recovers 60%–80% of halo descendants from their original progenitors over the redshift range 3 < z < 7, a sample purity that is not possible by other observational selection methods.

The purity of the sample can also verified by considering dark matter halo evolution under extended Press–Schechter modeling (e.g., Lacey & Cole 1993). This is illustrated in Figure 4. For an halo of mass M₁ = 10¹¹ M_☉ at z₁ = 0 (with number density n = 3 × 10⁻⁴ Mpc⁻³), we compute, following Trenti et al. (2008), the evolution of its descendant mass M₂ > M₁ at z₂ < z₁ (solid black line shows the median of the distribution, while the blue and green lines are 1σ and 2σ confidence intervals, respectively). The red squares in the figure show the mass at z₂ of dark matter halos at fixed number density n = 3 × 10⁻⁴ Mpc⁻³, demonstrating that number density-selected samples indeed trace the evolution of dark matter halos very well down to z ∼ 2, and extending the validity of the Papovich et al. (2011) findings to a wider redshift range. Interestingly, the number density selection appears to be less effective at z ≲ 2, although it still captures the growth of dark matter halos within 1σ.

Zoom In Zoom Out Reset image size

In order to create samples of constant cumulative number density in the HUDF, we begin by measuring the cumulative mass function of galaxies covering a wide range of redshifts and stellar masses, as shown in Figure 5. A compilation of these measurements is additionally provided in Table 2. Because of the exceptionally small size of the HUDF, cosmic variance contributes significantly to the error budget in this measurement. Errors have therefore been calculated according to the prescriptions of Moster et al. (2011), to account for the fluctuations in the observed galaxy number counts that are expected as a function of mass for a field of this size. We fit an exponential curve to each slice in redshift, which models the data well. The overlap of the curves tracing 2 < z ⩽ 3 and 3 < z ⩽ 4 in the high-mass regime (log M_*[M_☉] > 10.5) is likely one result of cosmic variance in this exceptionally small field.

Zoom In Zoom Out Reset image size

In Figure 6, we compare our measurements of the evolving number density of galaxies with logM_*[M_☉]  > 10 with published values from the literature. We find excellent agreement with the results of Marchesini et al. (2009) and González et al. (2011), whose measurements overlap in redshift with the available range of our sample. While the measurements of Marchesini et al. (2009) have been calculated assuming a Kroupa IMF (unlike those of González et al. (2011), which assume a Salpeter IMF), the errors have been calculated so as to accommodate variance due to this difference. Thus, these measurements are directly comparable and are fully consistent within the quoted errors. The literature data presented in this plot have likewise been corrected to a Salpeter IMF and compared using the most commonly applied cosmology (flat, Λ-dominated CDM with h₀ = 0.7 and Ω_m = 0.3).

Zoom In Zoom Out Reset image size

We define samples with constant cumulative number density (hereafter, Φ_c) for each redshift range by determining the stellar mass at the intersection of the exponential fits in Figure 5 and each of the horizontal dotted lines, equivalent to Φ_c = −3.0, −3.25, −3.5, and −3.75 Mpc⁻³. This provides the typical masses of galaxies and their progenitors at a constant number density over the range of observable redshifts. For each of the four number density-selected samples, we fit a linear function to the evolving log M_*[M_☉] at constant number density, as a function of redshift. We then select all galaxies in the photometric catalog with estimated stellar masses within 0.2 dex of the best-fit linear function at any given redshift.

For each of the chosen levels of Φ_c, in order of increasing number density, this method selects samples with total sizes of 14, 21, 40, and 53 galaxies, which will be later split into smaller bins by redshift for further analysis of their redshift evolution. Despite the unsurpassed depth of the WFC3/IR imaging, the small area of the HUDF field limits the number of samples with constant Φ_c that are fully accessible across the redshift range of 1 < z < 6. However, such measurements are possible for the galaxies with log Φ_c ⩽ −3.5 Mpc⁻³.

From an examination of the measured S/N of the IRAC and H₁₆₀ fluxes as a function of galaxy redshift and estimated stellar mass, we find these data are approximately complete for log M_*[M_☉] > 9.5 at z = 5. Thus, although the catalog includes galaxies at redshifts extending to z ∼ 8, none of our number density-selected samples are complete beyond z ∼ 6. Galaxies selected at log Φ_c ⩽ −3.5 are complete to a limit of z ∼ 6, but examining galaxies at progressively higher cumulative number densities requires a progressively lower mass detection threshold. Thus, our mass limit translates into progressively lower redshift limits for galaxies with increasingly higher cumulative number densities: galaxies selected at log Φ_c = −3.25 Mpc⁻³ are complete to z ∼ 5 and those at log Φ_c = −3.0 Mpc⁻³ are complete only to z ∼ 4.

Due to the exceptionally small area covered by the HUDF, our measurements become increasingly susceptible to cosmic variance for the most massive and lowest redshift galaxies. To circumvent this problem at low redshift, we have incorporated data from the K_s-selected FIREWORKS catalog of Wuyts et al. (2008) in order to provide a continuous measurement of the evolution of these number density-selected galaxy samples down to z ∼ 0.5. Although limited to lower redshifts, the FIREWORKS catalog spans a significantly larger region of overlapping sky, equal to 138 arcmin², and its included SEDs span an even broader range in wavelength coverage, extending from the U band to 24 μm MIPS imaging.

We calculate the photometric redshifts and stellar population estimates for galaxies in the FIREWORKS catalog with ⩾30% of maximum coverage in all observed bands, using the same EAZY and FAST prescriptions described in the previous sections. For objects included in both the H₁₆₀-selected catalog and FIREWORKS, comparisons of the photometric redshifts we determine with those provided by Wuyts et al. (2008) are shown in Figure 3. The photometric redshifts show close agreement, with little evidence for catastrophic failures and σ_NMAD = 0.055.

In Figure 7, we present the best-fit stellar masses of the individual H₁₆₀-selected galaxies as a function of redshift for four samples selected at different constant cumulative number densities. Linear fits to the moving average in each of the four samples are overplotted where the data are complete. We note that the data closely follow a linear relation in these plots as a direct result of the method of selection described in the previous section.

Zoom In Zoom Out Reset image size

The amplitude of each fit to the mass evolution at constant number density in Figure 7 rises with decreasing number density, as expected with this method of selection, indicating that samples at an increasing number density trace galaxies with decreasing mean stellar mass. The best fits to the slopes of the mass evolution do not significantly vary over the range in number density we examine, implying that the mass evolution occurs with approximately the same rate for each of these number density-selected samples. The result is an exponential build up equivalent to ∼1.5 orders of magnitude in stellar mass from z ∼ 6 to z ∼ 2 in each sample. These observations are consistent with well-established evidence showing that massive galaxies build up the bulk of their stellar mass prior to z = 2 (e.g., Thomas et al. 2005). The rate of growth we measure in the range 2 < z < 6 also precisely agrees with predictions for the mass evolution of star-forming galaxies from simulations (Oppenheimer & Davé 2008).

In Figure 8, we compare our results with the measurements of Papovich et al. (2011), which incorporated two separate methods for inferring the mass growth of galaxies with log Φ_c = −3.7 Mpc⁻³. The star-shaped points indicate measurements obtained using the integrated mass function of Marchesini et al. (2009) to determine the limiting mass of a sample with their desired cumulative number density. The square points indicate measurements of Papovich et al. (2011) that were indirectly determined using the integrated combined UV luminosity functions of Reddy & Steidel (2009), Bouwens et al. (2007), and Bouwens et al. (2011). Each of these measurements has been corrected to match estimates from a Salpeter IMF (+0.2 dex).

Zoom In Zoom Out Reset image size

The slope of the mass growth we measure with time agrees with the mass-selected measurements of Papovich et al. (2011) at the 1σ level. However, the amplitude of the fit we recover at a similar number density (log Φ_c = −3.75 Mpc⁻³) is consistently a factor of ∼2 higher in stellar mass than the UV-selected measurements over the overlapping redshift range. The method of rest-frame UV selection employed by Papovich et al. (2011) in the higher redshift regime may account for this difference, since such a selection is likely to be less inclusive of galaxies with lower SFRs.

Recent simulations by Krumholz & Dekel (2012) suggest that the mean mass estimates of UV-selected galaxies from Papovich et al. (2011) may suffer from >75% incompleteness in their sampling of star-forming central galaxies in the high-redshift regime. This incompleteness is suspected to be due to suppressed star formation scattering galaxies below the detection threshold at any given mass. As 45% of K-selected galaxies have been found to exhibit suppressed star formation at z = 2 (Kriek et al. 2006), Krumholz & Dekel (2012) predict a ∼50% incompleteness in mass-selected samples at high redshift.

For comparison with our measurements, we overplot the theoretical mass growth for a fully complete sample from Krumholz & Dekel (2012) as a solid cyan line in Figure 8. The two dashed lines with progressively lower amplitudes give the expected growth curves of Krumholz & Dekel (2012) for samples with 50% and 75% incompleteness, respectively. The mass evolution for a 50% complete number density-selected galaxy sample estimated by Krumholz & Dekel (2012) agrees with our observations to within 1σ. Thus, star formation duty cycle effects could yet be contributing to some degree of incompleteness in our sample, despite the apparent improvement over the measurements of Papovich et al. (2011) that we achieve by detecting in H₁₆₀ and incorporating SED constraints at 3.6 μm and 4.5 μm. The fact that our data exhibit a closer agreement with both observational analyses, as compared with the predictions of Krumholz & Dekel (2012), may also imply that the models require additional complexity to reproduce the results from observations.

The cosmological hydro-simulations of Jaacks et al. (2012) have also suggested that the duty cycle of bursty star formation histories may result in incompleteness in high-redshift observations. However, their simulations imply that this effect only significantly reduces the number of galaxies above the detection threshold of the HUDF for stellar masses below log M_*[M_☉] = 9.0. As the mass-completeness threshold of our analysis (log M_*[M_☉] > 9.5) is above this lower limit, it is not clear that bursty star formation histories alone can account for the difference in the number of galaxies observed, relative to the model expectations from Krumholz & Dekel (2012).

The evolution of the stellar masses selected at constant cumulative number density appears fairly smooth, largely as a consequence of the method of selection. However, the star formation properties of the same galaxies exhibit greater scatter, relative to that of their best-fit stellar masses, when selected by the same method. Thus, in order to examine the typical evolution in SFRs with better statistical certainty, in the following section we further bin the data by redshift to examine the mean evolution in SFR among the four sub-samples previously chosen.

In Figure 9, we present the redshift evolution of the averaged best-fit template SEDs for galaxies in our sample, selected to have log Φ_c = −3.5 Mpc⁻³. Each composite spectrum contains between 2 and 6 galaxies, the best-fit template SEDs of which have been shifted into the galaxy rest-frame and normalized at 0.75 μm. The evolution of the composite spectral shape over the redshift range 1 < z < 6 is consistent with an aging stellar population, providing a qualitative proof of concept for the method of cumulative number density selection and the usefulness of examining the mean evolution of stellar populations in galaxies selected by this method.

Zoom In Zoom Out Reset image size

In Figure 10, we present results for the cosmologically averaged evolution of SFRs in bins of δz = 1.25 for the four samples of galaxies, defined by cumulative number densities of log Φ_c = −3.0, −3.25, −3.5, and −3.75 Mpc⁻³. At redshifts where our galaxy samples are incomplete, we plot the mean SFR with empty circles with errors to indicate that these measurements should be regarded as upper limits. Incorporating measurements from the FIREWORKS catalog enables higher resolution of the peak SFR in each sample as well as the SFR evolution at z < 2.

Zoom In Zoom Out Reset image size

For each of the galaxy samples, we observe a steady increase in SFR from z ∼ 6, which turns over around 2 ≲ z ≲ 3. This shape is well approximated by an asymmetric Lorentzian function, defined by

$$\begin{equation} y=\frac{p_{0}p_{1}}{2}\times \frac{\frac{p_{1}}{2}+2p_{2}(x-p_{3})}{(x-p_{3})^{2}+(\frac{p_{0}}{2})^{2}}. \end{equation} \tag{ 5 }$$

We fit each of the four samples with this profile and shade the regions defined by a 68% confidence intervals in Figure 10. We stress that this shape has no obvious physical interpretation, particularly given the logarithmic nature of the plot. However, the χ² values returned by fitting the data to this form are approximately unity, indicating that this profile may be an acceptable function for interpolating through the points of measurement. While previous studies have chosen to fit a power law to the rising SFR over this redshift range, we find that the asymmetric Lorentzian function better models the more gradual rise of the data, as well as providing a smooth transition to the downturn below z ∼ 2.

The fits to the data exhibit no significant difference in the slopes of the SFR evolution in the range 2 < z < 6 among each of the samples selected at different constant number densities. Our fits also suggest that all four samples achieve a similar maximum SFR (log SFR [M_☉ yr⁻¹] ∼2) at 2 ≲ z ≲ 3. We find that galaxies with log Φ_c = −3.75 Mpc⁻³ reach a peak mean SFR at z ∼ 3; galaxies with −3.5 ⩽ log Φ_c ⩽ −3.25 Mpc⁻³ peak at z ∼ 2.3; galaxies with log Φ_c = −3.0 Mpc⁻³ peak at z ∼ 2. These results indicate that the timing of the peak SFR is correlated with the cumulative number density, such that samples with lower cumulative number density (and thus, higher average stellar mass) reach a peak SFR at progressively earlier times. This behavior appears consistent with what one would expect from cosmic downsizing (Cowie et al. 1996), in which the epoch of major mass build up in galaxies moves to lower redshift for increasingly lower-mass galaxies.

For comparison, we have overplotted measurements of the SFR evolution from Papovich et al. (2011) in Figure 10. Although the mean SFR estimates from Papovich et al. (2011) are approximately 2 times higher than what we measure at z > 3, they still agree within the stated 1σ errors. The evolution we measure in the SFR from 2 < z < 6 is thus consistent with the results of Papovich et al. (2011), despite the slight inconsistencies we find in the stellar mass evolution for the same cumulative number density (see also the comparison in Smit et al. 2012). This agreement in the measured SFR evolution is unsurprising, given that our rest-frame UV measurements are derived from overlapping data sets.

This rise in SFR with time, which we observe in each of the number density-selected samples prior to z ∼ 3, is predicted by theoretical models and hydro-simulations (e.g., Robertson et al. 2004; De Lucia et al. 2006; Finlator et al. 2006, 2007; Brooks et al. 2009; Trenti et al. 2010; Finlator et al. 2011; Krumholz & Dekel 2012). The number density-dependent evolution we measure in the cosmically averaged SFRs is also consistent with the recent simulations of Krumholz & Dekel (2012), who incorporate the effects of an evolving metallicity in modeling the star formation histories of galaxies. In Figure 10, we overplot the predictions of Krumholz & Dekel (2012) for galaxies selected over a similar range in cumulative number density. The agreement is impressive, especially for galaxies with log Φ_c ⩾ −3.5 Mpc⁻³, where our data are well sampled. The fact that the models closely track even the measurements at z > 5 suggests that the mean SFRs we measure in the extreme high-redshift regime are close to the predicted cosmological average, despite our possible mass incompleteness at those redshifts due to limitations from our photometric wavelength coverage and IRAC imaging depth. It is difficult to determine whether the departure of our measurements from the models of Krumholz & Dekel (2012) in the lowest number density sample (log Φ_c = −3.75 Mpc⁻³) is physical, as it could also be explained by the small number of galaxies included in the measurement, which will be most dramatically affected by cosmic variance.

With the stellar mass and SFR evolution each separately measured for the samples we have selected, we now consider how well the measured SFRs alone can account for the stellar mass build up we observe. In Figure 11, we use the best fits to the evolution of the cosmologically averaged SFRs to predict the stellar mass growth for the same samples of galaxies. In each panel, three tracks of mass growth derived from the integrated best-fit SFR evolution with respective duty cycles of 25%, 50%, and 100% are plotted alongside the separately measured mean stellar mass as a function of redshift. For each mass growth track with a duty cycle of 100%, we also provide the 68% confidence interval. In each case, we assume a fractional mass loss of 30% from star formation-driven winds, as traditionally applied to stellar populations modeled with a Salpeter IMF. We note that the true fractional mass loss during star formation will depend on a number of factors, such as the choice of initial mass function and recycling timescale.

Zoom In Zoom Out Reset image size

At all number densities we examine, the mass growth inferred from the SFR evolution broadly agrees with the observed mass evolution at z ∼ 2. However, at lower redshifts, we find that the highest mass galaxies in our sample appear to grow more massive than what is suggested from their observed SFRs alone. For galaxies at the lowest number density we examine, which are equivalent to galaxies with log(M[M_*]) ∼11.5 at z = 0, we find that even a duty cycle of 100% underpredicts the z = 1.5 stellar masses by more than a factor of three with >68% confidence. The inability of star formation alone to produce the observed build up in stellar mass in this high-mass sample indicates that additional processes are responsible for contributing the residual mass growth, particularly at late times (z < 2). Mergers are the most likely mechanism for making up this difference, although obscured star formation could also play a substantial role.

Among the galaxies selected at progressively higher number densities, we find an increasingly better agreement between the mass growth expected from the evolution of the observed SFRs within a reasonable range of duty cycles. The stellar mass build up observed for the samples with log Φ_c = −3.0 Mpc⁻³ can be fully accounted for until z ∼ 1.5 with a duty cycle as low as 25%. Galaxies with cumulative number densities in the range −3.5 ⩽ log Φ_c ⩽ −3.25 Mpc⁻³ are also found to agree with mass build up from star formation alone until z ∼ 2, assuming duty cycles in the range 50% −100%.

In general, at lower cumulative number densities (and thus, higher mean stellar masses), the mass growth from star formation alone falls increasingly further short of the observed stellar mass measurements. These data indicate that the fraction of mass build up that may be accounted for from star formation with a constant duty cycle scales with cumulative number density, such that galaxies with higher number density build up more of their mass through star formation alone, compared with their counterparts at lower number densities. Considering that the mean stellar mass of a galaxy is inversely related to its cumulative number density, it may not be surprising that the mechanism of stellar mass build up is dependent on the number density. Lower mass galaxies will have smaller gravitational potentials from which to attract mergers and recapture ejected mass from feedback processes. The fact that the gap between the observed stellar mass and that which is expected from the observed SFRs becomes especially prominent at z < 2 suggests that mergers contribute significantly to the late-time mass growth in the galaxies with log Φ_c ⩽ −3.5 Mpc⁻³. This finding is consistent with van Dokkum et al. (2010), who determined that in situ star formation can only account for a small fraction (∼20%) of the mass growth of galaxies with similar mean mass at z < 2.

At face value, our results may appear to conflict with the findings of Papovich et al. (2011), in which the observed SFRs were found to fully account for the average mass build up of galaxies with even slightly lower number densities than those examined here (log Φ_c = −3.7 Mpc⁻³). We note that the gap we find between the observed stellar mass growth and the growth inferred from the SFR only becomes significant at z < 2, the upper limit of the primary redshift range examined by Papovich et al. (2011). However, the fact that the Papovich et al. (2011) SFR evolution measurements were fit with an exponential function, which continually increases with time, suggests that the difference in fitting models could intrinsically produce greater build up of mass compared with ours in the low-redshift regime (z < 3), where the SFRs actually appear to decline.

In order to examine the effect of the chosen SFR evolution fitting function on the implied stellar mass growth, in Figure 11 we directly compare the stellar mass growth inferred from the Papovich et al. (2011) SFR evolution and stellar mass estimates, after fitting with the same functional form applied to our data and assuming a 100% duty cycle. In agreement with the results of Papovich et al. (2011), we find that the observed SFR evolution can fully account for the stellar mass growth, as estimated by Papovich et al. (2011) to z ∼ 3. At lower redshifts, their SFR evolution data become incomplete, but the trend of inferred mass growth suggests that the observed stellar mass growth begins to exceed the implied growth from the SFR alone z < 2, in agreement with our results. We refer the reader to Section 5.1 for a discussion of the differences in the observed stellar mass estimates.

In the previous two sections, we separately examined the cosmically averaged evolution of the mass and SFRs of cumulative number density-selected galaxies, providing new information about how galaxies grow, on average, in stellar mass over time. In the following section, we examine the relationship between SFR and stellar mass on the basis of individual galaxies, which provides a separate, valuable metric for understanding galaxy evolution.

Recent studies have shown that the SFRs of galaxies are inversely correlated with stellar mass at z < 2 (Zheng et al. 2007; Damen et al. 2009), consistent with models of cosmic downsizing. Observations at higher redshift have thus far indicated that sSFR of galaxies at a constant mass remains approximately constant with time from z ∼ 7 to z ∼ 2 (Stark et al. 2009; González et al. 2010; Labbé et al. 2010a, 2010b). These observations contrast the steady decrease in sSFR from high redshift, which is expected from the majority of theoretical models (Bouché et al. 2010; Davé et al. 2010; Dutton et al. 2010; Weinmann et al. 2011). Simulations of an evolving metallicity-dependent SFR have been evoked to explain the observed sSFR plateau at z > 2 (e.g., Krumholz & Dekel 2012). Improvements in the dust corrections of high-redshift observations have also been shown to be capable of raising the measured sSFR of galaxies z > 4 toward a closer agreement with theory (Bouwens et al. 2011). However, observations of the globally averaged sSFRs at z > 2 still continue to fall significantly short of the predictions from models.

It is worthwhile to determine whether the sSFR plateau at high redshift persists for galaxies selected at constant number density, rather than at constant mass. We thus present bin measurements of the sSFRs for the same four samples of galaxies selected at constant cumulative number density, shown as a function of redshift, in Figure 12. In agreement with the results of Noeske et al. (2007) and Daddi et al. (2007), we find that the globally averaged sSFR declines dramatically with time from z ∼ 2. Our results also indicate that this declining sSFR extends to z = 6, at least for the galaxies in the two lowest cumulative number density (and thus, highest mean mass) bins. Each of the number density-selected samples declines in sSFR from z ∼ 4 to z ∼ 2 and our data suggest that the trend may extend to z ∼ 8. However, mass constraints become increasingly difficult at such high redshifts, given our available wavelength coverage. Thus, the likely incompleteness in our data at z > 6 prohibits drawing definite conclusions from the apparent trend at those extreme redshifts.

Zoom In Zoom Out Reset image size

Intriguingly, the measurements for the two lowest number density samples indicate approximately one order of magnitude increase in the globally averaged sSFR from z = 2 to z = 6, precisely in agreement with the change in the specific gas accretion rate, which is predicted to scale as (1+z)^2.5 (Neistein & Dekel 2008). We have overplotted this predicted specific gas accretion rate, normalized to fit the amplitude of our measurements at z > 2 for each sample, in Figure 12. While incompleteness in our high-redshift measurements complicate the comparison for the samples with log Φ_c ⩾ −3.5 Mpc⁻³, the decline in sSFR agrees exceptionally well for the two lowest number density bins. In all cases, the observed mean sSFR begins to decline more steeply than the specific gas accretion rate model at times later than the peak epoch of star formation (z ≲ 2; see Figure 10).

In addition to the clear evidence for a decreasing sSFR from z ∼ 6 to z ∼ 2, we find a number density-dependent trend: at any given redshift, the sSFR is higher for galaxies at lower number density. This effect can qualitatively explain the fact that previous results (e.g., Stark et al. 2009; González et al. 2010; Labbé et al. 2010a, 2010b; Bouwens et al. 2011), which considered the evolution of the sSFR of galaxies at constant stellar mass, observe a plateau at z > 4, while we find a constantly declining sSFR from z ∼ 6. A constant mass selection will effectively probe lower number density samples with increasing redshift. As our findings suggest that the sSFR decreases with decreasing number density at all redshifts examined in this work, a plateau in the sSFR evolution would naturally result from a selection at constant mass. Thus, our findings are not in disagreement with previous results. Instead, we propose that the method of galaxy selection we employ better illuminates the sSFR evolution of matched galaxy samples over the same redshift range.

We note that our results are also consistent with the previous measurements of sSFRs measured for galaxies grouped cumulatively in mass at various redshift intervals in Damen et al. (2009). Still, due to large differences in the redshift and mass limits of the two samples, only one data point in our respective analyses may be directly compared. Damen et al. (2009) measure an average sSFR of 2 × 10⁻¹⁰ yr⁻¹ for galaxies with log(M_*/M_☉) > 11 at z = 1.5, precisely in agreement with the sSFR we measure for galaxies matched to the same mass and redshift interval.