THE SPECTRAL ENERGY DISTRIBUTIONS OF z ∼ 8 GALAXIES FROM THE IRAC ULTRA DEEP FIELDS: EMISSION LINES, STELLAR MASSES, AND SPECIFIC STAR FORMATION RATES AT 650 MYR^*

Deep imaging of high redshift z > 6 galaxies at mid-infrared (IR) wavelengths with the InfraRed Array Camera (IRAC; Fazio et al. 2004) on Spitzer has led to the surprising discovery that z ⩾ 6 galaxies have red H − [3.6] ∼ 0.5 colors. This has been taken as evidence of substantial stellar masses ∼10⁹–10¹⁰ M_☉ (Eyles et al. 2005; Yan et al. 2006; Stark et al. 2009) and stellar ages (>300 Myr; Labbé et al. 2006, 2010a, 2010b; González et al. 2010). Early studies also suggested that the specific star formation rate (sSFR) at fixed stellar mass was nearly flat at 3 < z < 8, in apparent disagreement with the strongly increasing specific inflow rates of baryons predicted by galaxy formation models (e.g., Neistein & Dekel 2008; Davé et al. 2011; Weinmann et al. 2011).

A key uncertainty in the analysis of the photometry is how much rest-frame optical emission lines contribute to the broadband fluxes. If emission lines are very strong, then the stellar masses previously derived by fitting stellar population models without lines would be biased high, and the sSFR would be biased low (e.g., Labbé et al. 2010b; Schaerer & de Barros 2010). To address this issue, emission lines have been included in some stellar population models (e.g., Schaerer et al. 2012; Schaerer & de Barros 2010; de Barros et al. 2012; Yan et al. 2012). As model predictions for emission lines are very uncertain empirical approaches are complementary and necessary.

Fumagalli et al. (2012) studied the evolution of the Hα equivalent width (EW) at 0 < z < 2 using near-IR spectra and predict very large values at z ∼ 8. If the emission lines are indeed as strong as argued, then deep photometry should be able to reveal the effect of the lines moving through the IRAC passbands with redshift. This technique has been applied successfully by looking for excess flux in the [3.6]-band due to Hα at z ∼ 4.5 (e.g., Shim et al. 2011, Stark et al. 2013) and by comparing the average observed IRAC colors over 4 < z < 6 (González et al. 2012a, 2012b). Nevertheless, the current situation is that the contribution of emission lines at 4 < z < 8 is still poorly known.

Direct estimates at z ∼ 8 would provide the best constraints yet on the evolution of emission line strengths and sSFR to the highest redshifts. Fortunately, at redshift z ∼ 8 a relatively clean test is possible for the strength of the strongest lines [O iii] λλ4959, 5007+Hβ, because these lie isolated in the [4.5] filter. Then main challenge is the extreme faintness of the sources.

In this Letter, we use the largest sample of z ∼ 8 galaxy candidates in combination with newly acquired ultradeep IRAC data from the IRAC Ultra Deep Field (IUDF) program (PI: Labbé) to study their colors, spectral energy distributions (SEDs), the contribution of emission lines, and to derive emission line corrected stellar masses and sSFRs. Throughout this paper, we assume an $\Omega _M=0.3,\Omega _\Lambda =0.7$ cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹. Magnitudes are in the AB photometric system (Oke & Gunn 1983).

The data analyzed here consist of ultradeep WFC3/IR imaging from the HUDF09 program (GO 11563; PI: Illingworth) over the Hubble Ultra Deep Field (HUDF) and two nearby fields HUDF09-1 and HUDF09-2, supplemented with deep WFC3/IR data observations from the Early Release Science program (GO 11359; PI: O'Connell) and the Multi-Cycle Treasury program CANDELS (PI: Faber/Ferguson; Grogin et al. 2011; Koekemoer et al. 2011) over the GOODS-South.

We use new ultradeep Spitzer/IRAC imaging from the IUDF program (PI: Labbé, PID 70145), a 262 hr Spitzer warm mission program at [3.6] and [4.5] μm (I. Labbé et al., in preparation). This survey increases the exposure time over the HUDF, HUDF09-1 and HUDF09-2 fields from 12–46 hr to ∼120 hr, ∼50 hr, and 80–120 hr, respectively. For the wider GOODS area we use the 23–46 hr deep IRAC coverage of GOODS (M. Dickinson et al., in preparation).

Our primary sample consists of 60 Y-dropout galaxies at z ∼ 8 selected by Bouwens et al. (2011) over the HUDF09 and Early Release Science (ERS) fields and 16 brighter Y-dropouts selected by Oesch et al. (2012) over the CANDELS GOODS-South area.⁸

We derive new IRAC photometry of all 76 sources following the procedure of Labbé et al. (2010a, 2010b). Briefly, we subtract nearby foreground sources based on their Hubble Space Telescope (HST) image profiles and determine local backgrounds (see, e.g., González et al. 2011). Then we perform aperture photometry through 2'' diameter apertures on the cleaned images and correct the flux for light outside the aperture using a profile derived from nearby stars (2.4 × in [3.6] and 2.5 × in [4.5]). Details of the photometry are presented in I. Labbé et al. (in preparation). We exclude 13 sources for which clean subtraction was not possible due to the proximity of very bright foreground sources, leaving a final sample of 63 Y-dropouts.

Figure 1 presents image stamps of IRAC-detected Y-dropout galaxies. A direct comparison to earlier GOODS observations demonstrates the clear improvement in sensitivity with the new IUDF data.

Zoom In Zoom Out Reset image size

We also include fainter, IRAC undetected, galaxies in the analysis by stacking 15'' × 15'' image stamps centered on the sources. On the stacks we subtract any residual background in a concentric annulus 5'' < r < 75. We compare between average and median combination to make sure that average results are not driven by outliers and derive uncertainties by bootstrap resampling.

To explore the observed properties of z ∼ 8 galaxies, we first present color–magnitude and color–color diagrams which contain diagnostic information about their stellar populations. Figure 2 (left panel) shows the H₁₆₀ − [3.6] versus H₁₆₀ color–magnitude diagram. The improved IRAC depth and inclusion of brighter candidates result in a direct detection with signal-to-noise ratio S/N([3.6] + [4.5]) ⩾ 3 of 20/63 sources (32%), of which 11 (17%) are detected at ⩾5σ.

Zoom In Zoom Out Reset image size

The average colors of IRAC detected (⩾3σ) galaxies are 〈H − [3.6]〉 = 0.52 ± 0.09, with an intrinsic scatter of σ(H − [3.6]) = 0.40 ± 0.12 (photometric uncertainties subtracted in quadrature), and 〈[3.6]–[4.5]〉 = 0.43 ± 0.11 with σ(H − [3.6]) = 0.55 ± 0.15. Stacking fainter H ∼ 27.3 (∼L*(z = 8)) galaxies, including IRAC undetected sources, shows similar colors. The mean 〈H − [3.6]〉 color is in the range of predictions from dust-reddened Bruzual & Charlot (2003, BC03) stellar population synthesis models, but the 〈[3.6]–[4.5]〉 appears significantly redder by ∼0.5 mag.

No obvious trend of H₁₆₀ − [3.6] color with H₁₆₀ or [3.6] mag is seen, but the H₁₆₀ − [3.6] (rest-frame U₁₇₀₀ − B₄₀₀₀ at z ∼ 8) does appear to correlate (at ∼2σ confidence) with J₁₂₅ − H₁₆₀ colors consistent with expectations from dust-reddened BC03 stellar population models. The trend is seen both in IRAC detected galaxies and in H₁₆₀ ⩽ 28 selected samples (⩽0.3L*(z = 8)), including IRAC undetected galaxies. In contrast, the observed [3.6]–[4.5] colors (rest-frame B₄₀₀₀ − V₅₀₀₀) are generally redder than the BC03 models, with larger offsets at bluer J₁₂₅ − H₁₆₀ colors.

Figure 3 shows the average SEDs of IRAC detected galaxies at z ∼ 8 compared those at z ∼ 7, only 130 Myr later. It is obvious that the observed [3.6]–[4.5] colors are markedly different. Shifted to the rest-frame, the SEDs viewed together suggest a clear flux excess at ∼5000 Å, shifting from [3.6] at z ∼ 7 to [4.5] at z ∼ 8. The excess is likely caused by [O iii]4959, 5007 and Hβ emission lines, which are expected to be strong in young star forming galaxies at low metallicities (e.g., Schaerer & de Barros 2010).

Zoom In Zoom Out Reset image size

Nevertheless, it is difficult to determine the emission line contribution using stellar population model fits to the broadband photometry at each redshift separately. Degeneracies between age and dust in the model colors and systematic uncertainties in implementation of emission lines hamper such a direct solution. Instead, we will attempt to circumvent this limitation by using the joint SED information at the two adjacent redshifts to isolate the effective emission line contribution.

Figure 4 shows clearly that the observed [3.6] − [4.5] colors at z ∼ 8 are 0.8 mag redder compared to z ∼ 7, while the H − [3.6] colors are 0.4 mag bluer.⁹ Changes in stellar population age and/or dust cannot produce such differences, but strong optical emission lines naturally reproduce the observed change.

Zoom In Zoom Out Reset image size

With the reasonable assumption that the average rest-frame spectrum is the same at both z ∼ 8 and z ∼ 7, we can solve directly for the rest-frame EW of these emission lines.

Considering only [O iii] λλ4959, 5007 and Hβ and applying the redshift selection functions of Oesch et al. (2012) and Bouwens et al. (2011), a linear fit to the colors produces W_{[O iii] λλ4959,5007+Hβ} = 500 Å. The redshift distribution accounts for the fact that emission lines may contribute to both [3.6] and [4.5]. Using the tabulated emission lines values of Anders & Fritze-v. Alvensleben (2003) for sub-solar metallicity 0.2Z_☉ (Erb et al. 2006; Maiolino et al. 2008), we also correct for the smaller integrated contributions of Hα, Hβ, Hγ, [N ii] and [O ii], leading to ${W}_{[{\rm O\,\scriptsize{III}}]\ {\lambda \lambda 4959,5007}+{\rm H}\beta }=670$ Å.

Note that this estimate is to first order independent of stellar population models and the fit only has one free parameter: the combined rest-frame EW of the strong emission lines. The main systematic uncertainty is the redshift selection function for the z ∼ 7 sample, as the [O iii] line moves from the [3.6] to the [4.5] filter at z > 7.0. A secondary uncertainty is the effect of metallicity, which changes the emission line ratios. To evaluate the impact of both effects we shift the central redshift of the z ∼ 7 distribution by up to Δz = 0.3 in either direction and vary the metallicity from 0.2–1.0 Z_☉. We add these systematic uncertainties in quadrature to the photometric errors arriving at ${W}_{[{\rm O\,\scriptsize{III}}]\ \lambda \lambda {4959,5007}+{\rm H}\beta }=670^{+260}_{-170}$ Å.¹⁰

The rest-frame EW implies a boost of $0.55^{+0.16}_{-0.11}$ mag to [4.5] and smaller $0.19^{+0.10}_{-0.09}$ mag to [3.6] at z ∼ 8, meaning the color of the stellar continuum is [3.6] − [4.5] ≈ 0. Using the line ratios of Anders & Fritze-v. Alvensleben (2003) for a metallicity of 0.2 Z_☉ and taking into account the slope of the continuum (f_ν ≈ constant or f_λ∝λ⁻²) the corresponding Hα equivalent is ${W}_{{\rm H}\alpha }=430^{+160}_{-110}$. Note that the inferred WHα depends on assuming line ratios.

A common approach to estimate the sSFR is derive it from W_Hα, using the Hα to star formation rate (SFR) conversion of Kennicutt (1998) and the stellar continuum from BC03 stellar populations synthesis models. The model adopts a Salpeter (1955) initial mass function between 0.1–100 M_☉, a metallicity of 0.2 Z_☉, and constant SFR. Using the indirect estimate of Hα we derive ${\rm sSFR}=11^{+11}_{-5}$ Gyr⁻¹, corresponding to an approximate stellar age of $100_{-50}^{+100}$ Myr. The estimated sSFR is not very sensitive to the assumed star formation history (e.g., declining and rising SFH introduce changes <0.1 dex), but potentially significant uncertainties are the intrinsic line ratios, the possibly larger attenuation of the emission lines compared to the continuum and the escape fraction of ionizing photons. Larger attenuation and a non-zero escape fraction would increase the derived sSFR.

Fitting BC03 models with CSF to the broadband SEDs using FAST (Kriek et al. 2009), we inspect the stellar masses and mass-to-light ratios after correcting the IRAC [3.6] and [4.5] bands for emission line contributions. The derived mass-to-light ratios of IRAC detected galaxies at z ∼ 8 are M/L₁₈₀₀ ≈ 0.05 and M/L_V = 0.09, in good agreement with those found by Yan et al. (2012) for two bright z ∼ 8 sources. The M/L would have been overestimated by ∼3 × if emission lines had not been corrected. The scatter in M/L, based on the intrinsic scatter in H − [3.6] and assuming CSF, is σ(M/L₁₈₀₀) = 0.25 dex and σ(M/L_V) = 0.1 dex, where the luminosities refer to the stellar continuum only. This is likely a lower bound to the true scatter in M/L due to selection effects. The typical L*(UV, z = 8) galaxy (corresponding to H ∼ 27.3) has a median rest-frame color U − V ∼ 0.2 (stellar continuum only) and a stellar mass $M=0.9_{-0.4}^{+0.6}\times 10^{9} \,M_\odot$.

Model fits to the corrected SEDs can also be used to determine sSFRs, although these are rather uncertain due to the degeneracies between age and dust and the weak constraints on the UV–slope β at z ∼ 8 (Bouwens et al. 2013). To mitigate the difficulties somewhat, we can fix the reddening to A_V = 0.40 ± 0.15, consistent with the average β values of bright z ∼ 7 galaxies (Bouwens et al. 2013). The best-fit sSFR is then ${\rm sSFR}=7_{-4}^{+13}$ Gyr⁻¹, in good agreement with the sSFR derived from W_Hα.

We derive integrated stellar mass densities of H-band selected galaxies to faint UV-limits by converting the stepwise UV–luminosity function (LF; Bouwens et al. 2011; Oesch et al. 2012) to stepwise mass densities, using the best-fit M/L to the emission line corrected SED.

Note that the emission line correction is derived from ∼1.5L*(UV, z = 8) galaxies, so we assume that the effect of emission lines is similar for L < L* galaxies. This approximation is unavoidable as sub-L* galaxies are currently beyond the reach of Spitzer, even at our ultradeep limits. We derive uncertainties by randomly drawing 10⁴ LFs, perturbing the stepwise LF and M/L ratios by their errors, and integrating the resulting mass densities. Integrating to M_{UV, AB} = −18 yields $\rho ^*(z=8)=0.6^{+0.4}_{-0.3}\times 10^6 \,M_\odot$ Mpc⁻³.

This paper presents new ultradeep IRAC data from the IUDF program (120 hr in [3.6] and [4.5] over the HUDF) in combination with a large sample of 63 Y₁₀₅-dropout galaxies at z ∼ 8 to study average HST/WFC3 and Spitzer/IRAC SEDs, the effect of emission lines, and their impact on stellar population model fits.

The most exciting result is the direct detection of large numbers of z ∼ 8 galaxies by Spitzer/IRAC to faint limits, without signs of a hard "confusion limit." The ability, for the first time, to determine the distribution of rest-frame optical colors at z ∼ 8 is a clear reminder of the enduring capabilities of Spitzer, currently in its post cryogenic ("warm") mission.

Using the joint IRAC constraints on the SED at z ∼ 7 and z ∼ 8, we derive an empirical constraint on the emission line contribution of rest-frame EW ${W}_{[{\rm O\,\scriptsize{III}}]4959,5007+{\rm H}\beta }\sim 670$ Å contributing ∼0.5 mag to [4.5] at z ∼ 8. To place this in context, we compare to recent estimates for star forming galaxies from other surveys (shown in Figure 5). While direct comparisons are difficult due to different techniques and stellar mass ranges, our results confirm the picture of an increasing importance of emission lines toward higher redshift (e.g., Shim et al. 2011; Fumagalli et al. 2012; Stark et al. 2013). The evolution of the rest-frame EW of Hα follows W_Hα∝(1 + z)^{1.2 ± 0.25} over 1 < z < 8. This a factor of ∼2 slower than the power law derived by Fumagalli et al. (2012) over 0 < z < 2, but in agreement with their model which assumes that the cumulative number density of progenitors at higher redshift remains constant (e.g., van Dokkum et al. 2010; Papovich et al. 2011).

Zoom In Zoom Out Reset image size

Clearly, at such high EW emission lines cannot be ignored when deriving stellar masses and ages at high redshift (e.g., Schaerer & de Barros 2010; de Barros et al. 2012). After correcting the SEDs for emission lines we find that the L*(z = 8) galaxies have stellar masses and M/L that are ∼3 × lower.

The W_Hα based sSFR$=11^{+11}_{-5}$ Gyr⁻¹ are ∼5 × higher than those at z = 2 (Reddy et al. 2012), consistent with an increase ${\rm sSFR}\propto (1+z)^{1.5^{+0.7}_{-0.5}}$ between 2 < z < 8. This evolution is somewhat higher than derived by González et al. (2012b) over 2 < z < 6 and somewhat lower than inferred by Stark et al. (2013) at 2 < z < 7, while in general agreement with the cumulative number density model of Fumagalli et al. (2012). The latest estimates reduce the tension between the observations and simulations, which predicted an evolution of the baryonic mass inflow rates ∝(1 + z)^2.3 (e.g., Neistein & Dekel 2008; Davé et al. 2011; Weinmann et al. 2011).

Obviously, a huge leap in the determination of the properties of galaxies at these redshifts will be made possible by the arrival of James Webb Space Telescope. But in the mean time wider and still deeper future IRAC observations will continue to provide valuable insight into the properties of the z ∼ 8 population, with emerging detections of z = 9–10 (e.g., Zheng et al. 2012) galaxies probably not far behind.

We are grateful to Jarle Brinchmann, Steve Finkelstein, Casey Papovich, Dan Stark, Simone Weinmann for helpful discussions. We thank the anonymous referee for useful suggestions which improved the paper. Support for this work was provided by NASA through Hubble Fellowship grant HF-51278.01. Support for this work was provided by NASA through an award issued by JPL/Caltech. This work has further been supported by NASA grant HST-GO-11563.01. We also acknowledge funding from ERC grant HIGHZ No. 227749.

Facilities: Spitzer (IRAC) - Spitzer Space Telescope satellite, HST (WFC3/IR) - Hubble Space Telescope satellite