The Brightest z ≳ 8 Galaxies over the COSMOS UltraVISTA Field

In an effort to further ascertain the nature of the z ∼ 8 LBG sample considered in this work, we also inspected the recent Drift And SHift mosaic (DASH—Momcheva et al. 2016; Mowla et al. 2019) at the nominal locations of the selected candidate bright LBGs. This mosaic covers ∼0.7 sq. deg of sky in the WFC3/F160W band to a depth of ∼25.1 mag (03 diameter aperture—Mowla et al. 2018), and overlaps approximately with three of the four UltraVISTA ultradeep stripes (see Figure 1). As a bonus, the mosaic also incorporates all the publicly available imaging in the F160W band over the COSMOS/UltraVISTA field. Given the detection of the candidate LBGs was performed on ground-based data (seeing FWHM ∼ 07), the finer spacial resolution of HST/WFC3 (PSF FWHM ∼ 02) is key to test potential multiple components of the candidate bright LBGs, whose blending could artificially increase their measured luminosity (e.g., Bowler et al. 2017; Marsan et al. 2019) or systematically affect their redshift estimates.

We found that 14 of the 16 candidate LBGs are covered by the DASH mosaic. Their image stamps are presented in Figure 3, while in Table 2 we summarize the coverage details for each source. We note that two sources (UVISTA-Y4 and UVISTA-Y8) fall on or very close to the border between the DASH coverage and deeper WFC3 coverage, resulting in unreliable measurements.

Zoom In Zoom Out Reset image size

Inspection of the DASH mosaic at the locations of the candidate LBGs discussed in this work resulted in single, isolated sources (for the five sources that are detected at ≳4σ), with the important exception of one candidate, UVISTA-Y3. In Figure 4 we present an image stamp extracted from DASH overplotted with the contour of the combined J, H, and K_s imaging data. A SExtractor run identified three individual objects (with S/N ∼ 4.5, 2.9, and 2.2) overlapping with the UltraVISTA footprint of UVISTA-Y3, that we label as UVISTA-Y3a, UVISTA-Y3b, and UVISTA-Y3c, for the three components in order of increasing decl., respectively (see Figure 4). The three sources are found to have relative distances of ∼05. To further ascertain the multiple nature of this source, we run a Monte Carlo simulation, presented in Appendix A, consisting in adding to the DASH footprint synthetic sources whose morphologies are similar to those measured for bright z ≳ 6 LBGs. None out of the 20 synthetic sources were split into multiple components by the background noise, increasing our confidence in the multi-component nature of this source. The high resolution provided by the DASH imaging enabled rerunning the photometry with mophongo, this time adopting the DASH image itself as positional and morphological prior. As we will show in the next section, the single z ∼ 8 source initially identified on the UltraVISTA images resulted in the three objects being at z ≳ 8.

Zoom In Zoom Out Reset image size

The relatively low S/N significance of the detections of the three components prevents from a comprehensive assessment of their morphology and associated uncertainties. A number of works have found that the typical effective radii for LBGs of luminosities similar to those in our sample and at similar redshifts are r_e ≲ 1 kpc (e.g., Holwerda et al. 2015; Oesch et al. 2016; Bowler et al. 2017; Stefanon et al. 2017b; Bridge et al. 2019). At z ∼ 8, a separation of 05 corresponds to ∼2.5 kpc, i.e., ≳2.5× the typical size of bright LBGs at these redshifts. In the spirit of providing further context, we performed an estimate of the sizes for the three sources using the method of Holwerda et al. (2015), and found effective radii of r_e ∼ 0.6, 0.5 and 0.4 kpc, respectively for UVISTA-Y3a, UVISTA-Y3b, and UVISTA-Y3c, further supporting our interpretation as three distinct sources. We stress though, that our r_e estimates are only indicative, and should not be considered out of this context.

In our deblending, the flux density of UVISTA-Y3b in the IRAC bands results to be marginal compared to that of the other two components. One possible explanation for this is that while UVISTA-Y3a and UVISTA-Y3c lie at opposite locations with respect to the observed peak of flux density, UVISTA-Y3b is offset from that. In such a configuration, the observed peak of flux density does not coincide with any of the detected sources; instead, it is likely the result of the overlap of the wings of the light profiles of these two components, suggesting the two sources could account for most of the observed flux density. To test this interpretation, we forced the exclusion of either UVISTA-Y3a or UVISTA-Y3c in the deblending process. The result was residual flux at the location of the corresponding component, suggesting these two sources are required to fully account for the observed IRAC flux. However, for a more robust determination of the deblended flux density, higher S/N observations with HST/WFC3 and possibly at wavelengths >3 μm (and/or higher spatial resolution) are likely needed. We therefore cannot be sure that our best-fit decomposition is entirely free from systematic errors.

Given that there are 16 z ∼ 8 candidates over the ∼0.8 deg² of the UltraVISTA ultradeep stripes, we would expect to find only ∼1 candidate over the ∼190 arcmin² CANDELS COSMOS field. Indeed, only one z ∼ 8 candidate from our selection is located over the CANDELS COSMOS field (UVISTA-Y11). In Figure 5 we present the image stamps in the V₆₀₆ I₈₁₄, J₁₂₅, JH₁₄₀, and H₁₆₀. The V₆₀₆ mosaic shows a close low-z neighbor just  ∼ 07 west of UVISTA-Y11, which is not detected in any NIR image (see Figures 5 and 6). Therefore, we manually included this low-z neighbor when performing the photometry.¹⁴ We do not detect flux at >1σ in the V₆₀₆ and I₈₁₄ bands, increasing our confidence in its high-z nature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we inspected the ACS I₈₁₄-band mosaic of the COSMOS program (Scoville et al. 2007, ∼26.5 mag in 06 aperture diameter, 5σ). We found coverage for all sources with the exception of UVISTA-Y1 and UVISTA-Y15. No significant detections exist for any of the sources. We identified a potential low-z galaxy ∼10 northwest of the nominal location of UVISTA-Y4, which, however, does not affect our flux density estimates.

The above analysis based on serendipitous deep HST coverage for two among the brightest z ∼ 8 LBGs stresses the importance of deep (≳1 orbit) high-resolution multi-band follow-up to further assess the nature of the remarkable LBG candidates identified in the present work.

Figure 6 presents the image stamps of all the candidate z ∼ 8 LBGs. Their positions and main photometry are listed in Table 3, while in Appendix B we list the flux densities for all objects in all bands. As it is evident from Figure 6, the majority of the sources are clearly detected in the near-infrared, and most of them are also detected in at least one of the Spitzer/IRAC bands. The brightest source has an H-band magnitude of 24.8 mag and it is detected at 12σ, adding in quadrature the detection significance in the J, H, and K_s bands.

Table 3.  Sample of Candidate z ∼ 8 LBGs

ID	R.A.	Decl.	mHa	Y − Jb	[3.6]–[4.5]b	zphotc
	(J2000)	(J2000)	(mag)	(mag)	(mag)
UVISTA-Y1d,i	09: 57: 47.900	+02: 20: 43.66	24.8 ± 0.1	>2.1	0.4 ± 0.2	${8.53}_{-0.62}^{+0.53}$
UVISTA-Y2i	10: 02: 12.558	+02: 30: 45.71	24.8 ± 0.2	>2.2	0.5 ± 0.1	${8.21}_{-0.49}^{+0.50}$
UVISTA-Y3ae	10: 00: 32.324	+01: 44: 30.86	25.5 ± 0.3	>0.9	0.6 ± 0.5	${8.68}_{-1.21}^{+0.93}$
UVISTA-Y3be	10: 00: 32.317	+01: 44: 31.48	26.1 ± 0.5	>0.9	<0.8f	${8.90}_{-1.18}^{+1.24}$
UVISTA-Y3ce	10: 00: 32.350	+01: 44: 31.73	26.0 ± 0.5	>−0.5	0.7 ± 0.5	${9.29}_{-2.10}^{+1.58}$
UVISTA-Y4i	10: 00: 58.485	+01: 49: 55.96	24.9 ± 0.2	1.0 ± 0.4	0.1 ± 0.2	${7.42}_{-0.20}^{+0.19}$
UVISTA-Y5d,i	10: 00: 31.886	+01: 57: 50.23	24.9 ± 0.2	>1.3	0.8 ± 0.3	${8.60}_{-0.65}^{+0.58}$
UVISTA-Y6d	10: 00: 12.506	+02: 03: 00.50	25.3 ± 0.3	>1.5	0.3 ± 0.4	${8.32}_{-0.92}^{+0.66}$
UVISTA-Y7i	09: 59: 02.566	+02: 38: 06.05	25.5 ± 0.4	>1.3	⋯h	${8.47}_{-0.73}^{+0.72}$
UVISTA-Y8i	10: 00: 47.544	+02: 34: 04.84	25.4 ± 0.3	>1.4	1.0 ± 0.8	${8.34}_{-0.58}^{+0.60}$
UVISTA-Y9	09: 59: 09.621	+02: 45: 09.68	25.4 ± 0.3	0.8 ± 0.7	⋯h	${7.69}_{-0.71}^{+0.99}$
UVISTA-Y10i	10: 01: 47.495	+02: 10: 15.37	25.3 ± 0.3	>1.6	0.9 ± 0.7	${8.25}_{-0.60}^{+0.61}$
UVISTA-Y11i	10: 00: 19.607	+02: 14: 13.15	25.2 ± 0.3	>1.4	0.8 ± 0.4	${8.64}_{-0.72}^{+0.66}$
UVISTA-Y12i	10: 00: 15.975	+02: 43: 32.96	25.6 ± 0.4	>1.2	0.2 ± 0.8f	${8.70}_{-0.74}^{+0.61}$
UVISTA-Y13	09: 58: 45.561	+01: 53: 41.79	25.8 ± 0.4	>1.1	0.8 ± 0.7	${8.54}_{-1.18}^{+0.79}$
UVISTA-Y14	10: 00: 12.568	+01: 54: 28.50	25.6 ± 0.4	>1.1	0.1 ± 0.6	${7.55}_{-2.68}^{+1.71}$
UVISTA-Y15	09: 57: 35.795	+02: 11: 57.81	25.6 ± 0.4	1.1 ± 0.9	<−0.5f,g	${7.64}_{-1.13}^{+1.13}$
UVISTA-Y16i	10: 01: 56.333	+02: 34: 16.25	25.3 ± 0.3	1.2 ± 0.7	0.6 ± 0.4	${7.90}_{-0.57}^{+0.74}$

Notes. Measurements for the ground-based bands are 12 aperture flux densities after removing neighboring sources with mophongo and corrected to total using the PSF and luminosity profile information; measurements for Spitzer/IRAC bands are based on 18 aperture flux densities from mophongo corrected to total using the PSF and luminosity profile information. We refer the reader to Tables 7–9 in Appendix B for a complete and more detailed listing of the flux density measurements for all objects in our sample.

^aH-band magnitude and associated 1σ uncertainty estimated from the UltraVISTA DR3 mosaic. ^bUpper/lower limits to be intended as 1σ. In computing these colors, we replaced negative fluxes with their corresponding 1σ uncertainty. See Tables 7–9 for a complete listing of flux densities in all bands. ^cPhotometric redshift and 68% confidence interval of the best-fitting template from EAzY. ^dThese sources were already presented in Stefanon et al. (2017b). We propose them here again for completeness, noting that their associated parameters in the present work were computed excluding the information from the HST bands. We refer the reader to Stefanon et al. (2017b) for a more complete analysis. ^eThese candidate LBGs were initially identified as a single source on the UltraVISTA NIR bands. Successive analysis including COSMOS/DASH suggests these are three distinct objects. The corresponding observables when a single object is assumed are: R.A. = 10:00:32.322; decl. = 1:44:31.26, m_H = 25.0 ± 0.1 mag; Y − J = 1.1 ± 0.4 mag; [3.6]−[4.5] = 0.3 ± 0.1 mag and ${z}_{\mathrm{phot}}={7.62}_{-0.28}^{+0.14}$. ^fThis IRAC color is based on <2σ flux density estimate in both bands. ^gA blue [3.6]−[4.5] < 0 mag color might be indicative of a redshift z ≲ 7. ^hAfter visual inspection, the neighbor-clean image stamps in the IRAC 3.6 and 4.5 μm bands, which constitute the base for our flux density estimates, showed non-negligible residuals that likely systematically affect our estimates. We therefore opted for excluding from our analysis the measurements involving IRAC for these sources. ⁱThese sources have a probability p(z > 7) ≥ 0.97, suggesting these may be fairly robust candidates of bright LBGs.

Download table as:  ASCIITypeset image

The observed SEDs of the galaxy candidates are presented in Figure 7, along with the EAzY best-fit templates at z ∼ 8 and, to provide contrast, forced fits to model z < 6 galaxies. The inset in each panel presents the redshift likelihood distribution based on the available optical, infrared, and Spitzer/IRAC photometry. Finally, in Figure 8 we show the SED of UVISTA-Y3 when we do not deblend its photometry using the information from the DASH imaging. This SED is best-fitted by a z ∼ 8 solution, consistent with our initial selection.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Four of our 16 z ∼ 8 candidates (or ∼23% of our sample) are located outside the region with the deepest optical observations from the CFHT legacy deep survey. Because the HSC imaging was not available at the time of the initial sample selection, and given shallower optical observations available in some of the bands to control for contamination (e.g., in the z-band), we can ask whether we find an excess of sources over these regions compared to what we would expect from simple Poissonian statistics. As the outer region contains ∼37% of the area, we find no evidence for a higher surface density of z ∼ 8 candidate galaxies outside those regions providing the best photometric constraints. This suggests that we can plausibly include the full UltraVISTA search area in quantifying the volume density of bright z ∼ 8 galaxies. Furthermore, the subsequent addition of flux densities from the HSC mosaics did not substantially affect the redshift distributions for these sources, increasing our confidence on their being at z ≳ 8.

Although most of our sample sources are robust z > 8 candidates, a few have relatively unconstrained redshift probability distributions. Specifically, 10 sources (when considering UVISTA-Y3 as multiple objects) have a 97% or higher probability of being genuine LBGs at z > 7, namely UVISTA-Y1, UVISTA-Y2, UVISTA-Y4, UVISTA-Y5, UVISTA-Y7, UVISTA-Y8, UVISTA-Y10, UVISTA-Y11, UVISTA-Y12, and UVISTA-Y16, while the remaining 8 sources, UVISTA-Y3a, UVISTA-Y3b, UVISTA-Y3c, UVISTA-Y6, UVISTA-Y9, UVISTA-Y13, UVISTA-Y14, and UVISTA-Y15, have probabilities 0.6 ≲ p(z > 7) ≲ 0.95.¹⁵ These tend to have the reddest J − H colors and hence the least certain breaks. Encouragingly enough, the most uncertain sources are distributed fairly uniformly across the UltraVISTA search area and are not located exclusively over those regions with the poorest observational constraints.

While 14 out of the 16 candidates do not present any significant detection in the 5.8 μm and 8.0 μm bands, two sources in our z ∼ 8 selection (UVISTA-Y3 and UVISTA-Y13) are formally detected at >1σ in the combined 5.8 μm and 8.0 μm observations, with nominal brightnesses of ∼23–23.5 mag at >5 μm. This could be interpreted as indication of contamination from intrinsically red z < 3 galaxies; however, assuming an intrinsic flux density of 350 nJy (∼25 mag, i.e., an approximately flat f_ν SED) at ∼7μm, simple noise statistics predict 4 ± 2 sources to be detected at >1σ. We therefore conclude that the >1σ formal detection of two z ∼ 8 candidates in our selection is not a concern.

The selection criteria expressed by Equations (2) and (3) are designed to select z ≳ 9 LBG candidates. Indeed our initial analysis identified two exceptionally bright (m_H ∼ 22.5 mag) J-dropouts (UVISTA-J1 and UVISTA-J2). However, follow-up analysis including our HST/WFC3 data and presented in Stefanon et al. (2017b) revealed that these two sources are likely z ∼ 2 interlopers. For this reason, we omit them from the present sample and refer the reader to Stefanon et al. (2017b) for full details.

In summary, to facilitate the comparison of our results to both simulations and observations of LBGs at z ∼ 8, in the rest of this work we consider the 16 Y-band dropouts as our fiducial sample of galaxies at z ∼ 8; specifically, we include in the z ∼ 8 sample those Y-dropouts with nominal z_phot ∼ 9 (see Figure 2). However, in Section 6.5 we also consider the contribution of those sources with z_phot ∼ 9 to the z ∼ 9 LF. We refer the reader to our discussion in Section 6.5 for full details.

One potentially important source of contamination for our current z ∼ 8 and z ∼ 9 samples occurs through the impact of noise on the photometry of foreground sources in our search fields. While noise typically only has a minor impact on the apparent redshift of various foreground sources, the rarity of bright z ∼ 8–10 galaxies makes it possible for the noise to cause some lower-redshift galaxies to resemble high-redshift galaxies similar to those we are trying to select. This issue tends to be most important for very wide-area surveys where there exist large numbers of sources that could scatter into our input catalog.

To determine the impact that noise can have on our samples, we started with an input catalog of z ≤ 6 sources (13000 in total) extracted from the CANDELS/3D-HST catalogs (Skelton et al. 2014; Momcheva et al. 2016) over the deep regions in the GOODS North and GOODS South fields, and with apparent magnitudes ranging from H₁₆₀ = 23 to 26 mag. The procedure was replicated 25 times randomly varying the flux densities according to the measured uncertainties to increase the statistical confidence and to simulate the expected number of sources in the 3000 arcmin² of the UltraVISTA field.

Fitting the photometry of each source to a redshift and the SED template set described in Section 4, we derived an SED model for each source in the catalog based on the available photometry and the EAzY SED templates. We then used that to estimate the equivalent flux for each source in the ground-based imaging bands available over UltraVISTA and perturbed those model fluxes according to the measured noise over the shallow and deep regions over UltraVISTA and according to the depth available over SPLASH, SEDS, and SMUVS. Finally, we reselected sources using the same selection criteria as we applied to the actual observations. In perturbing the fluxes of individual sources, we considered both Gaussian and non-Gaussian noise (the latter of which we implemented by increasing the size of noise perturbations by a factor of ∼1.3).

Our simulations suggested a very low contamination fraction for our z ∼ 8 samples. Over the ultradeep stripes where 95% of the sources in our z ∼ 8 sample were found, these simulations predicted just one z < 6 contaminant for the entire ∼0.8 sq. deg. area, equivalent to a contamination fraction of 5% for our z ∼ 8 samples. The typical H-band magnitude of the expected contaminants ranged from H ∼ 25 to 25.5 mag.

A number of recent works has shown that gravitational lensing from foreground galaxies could have a particularly significant effect in enhancing the surface density of bright z ≥ 6 galaxies (e.g., Wyithe et al. 2011; Barone-Nugent et al. 2015; Fialkov & Loeb 2015; Mason et al. 2015). This is especially true for the brightest sources due to the intrinsic rarity and the large path length available for lensing by foreground sources. It has thus become increasingly common to look for possible evidence of lensing amplification in samples of z ∼ 6–10 LBGs (e.g., Bowler et al. 2014, 2015; Oesch et al. 2014; Zitrin et al. 2015; Bernard et al. 2016; Bouwens et al. 2016; Roberts-Borsani et al. 2016; Morishita et al. 2018; Ono et al. 2018).

Even though the fraction of lensed sources among bright samples does not seem to be particularly high (Bowler et al. 2014, 2015), we explicitly considered whether individual sources in our bright z ∼ 8 galaxy compilation showed evidence for being gravitational lensed. For convenience, we used the Muzzin et al. (2013) catalogs providing stellar mass estimates for all sources over the UltraVISTA area we have searched. These catalogs use the diverse multi-wavelength data over UltraVISTA, including GALEX near- and far-ultraviolet, HST optical, near-infrared, Spitzer/IRAC, and ground-based observations, to provide flux measurements of a wide wavelength range and then use these flux measurements to estimate the redshifts and stellar masses. We also verified that the values obtained did not differ substantially (≲15%) from those obtained adopting the stellar mass estimates of Laigle et al. (2016).

As in Roberts-Borsani et al. (2016), we model the foreground objects as singular isothermal spheres (SIS) to assess their influence on the z ∼ 8 galaxy luminosities, and we use the measured half-light radius (Leauthaud et al. 2007) and inferred stellar mass to derive velocity dispersion estimates for individual galaxies in these samples. For cases where size measurements were not available from HST I₈₁₄-band imaging over the COSMOS field, we estimated the half-light radius relying on the mean relation derived by van der Wel et al. (2014). Of the 16 z ∼ 8 in our primary sample, only four appear likely to have their flux boosted (>0.1 mag) by lensing amplification.

One of the main advantages of the SIS model is the availability of analytic expressions for the main observables (e.g., magnification, shear, convergence) at the expense of a simplified (spherically symmetric) gravitational potential. For all of our candidate LBGs with the exception of Y6, the lenses have compact, quasi-spheroidal morphology (minor-to-major axis ratio b/a ≳ 0.9) supporting the adoption of a SIS model. For Y6, instead, of three lensing sources, only one has a spheroidal morphology, while the remaining two have elongated shapes (b/a ∼ 0.5), with a position angle of the LBG relative to the main axes of the two ellipses of ∼296 and ∼45, respectively.

More realistic magnification factors could be obtained for Y6 assuming a singular isothermal ellipsoid model (SIE—e.g., Kormann et al. 1994; Kochanek et al. 2004) for the two elongated lensing galaxies. In particular, if the major axis of the ellipsoid is oriented toward the high-redshift source, the magnification from a SIE model could be sensibly higher than the magnification from a SIS model. For the two elliptical lenses, the magnifications from the SIE model are 4.5% and 8.4% higher than the corresponding estimates from the SIS model, corresponding to ∼0.08 and ∼0.04 mag difference. Given the small contribution to the magnification estimates, and because the increase in magnification relative to SIS are just a fraction of the systematic uncertainties from the stellar mass estimates of the lensing sources (∼15%), in this work we adopt magnification factors from the SIS model for all lenses.

In the following, we present in more detail our estimates of lensing magnification for the four sources:

UVISTA-Y6: This source is estimated to be amplified by ∼1.4×, ∼1.16×, and ∼1.14× from a 10^10.7 M_⊙, z = 1.76 galaxy (10:00:12.51, 02:02:57.3), 10^10.6 M_⊙, z = 1.6 galaxy (10:00:12.15, 02:02:59.6), and a 10^10.3 M_⊙, z = 1.65 galaxy (10:00:12.18, 02:03:00.7), respectively, that lie within 49, 32, and 54 of this source. Their velocity dispersions are estimated to be 259 km s⁻¹, 225 km s⁻¹, and 206 km s⁻¹, respectively.

UVISTA-Y8: This source is estimated to be amplified by 1.39× from a 10^10.8 M_⊙ (264 km s⁻¹), z = 1.33 galaxy (10:00:47.68, 02:34:08.4) that lies within 41 of this source.

UVISTA-Y9: This source is estimated to be amplified by 1.37× and 1.43× by a 10^11.0 M_⊙ (265 km s⁻¹), z = 0.91 galaxy (09:59:09.35, 02:45:11.8) and 10^11.0 M_⊙ (268 km s⁻¹), z = 0.93 galaxy, respectively, that lie within 50 and 46 of the source.

UVISTA-Y13: This source is estimated to be amplified by 1.6× by a 10^11.15 M_⊙ (330 km s⁻¹), z = 1.63 galaxy (09:58:45.83, 01:53:40.6) that lies within 42 of the source.

We discuss the potential impact of lensing on our inferred value for the characteristic magnitude of the UV LF, M*, at the end of Section 6.7.

In Figure 9 we present our sample of candidate z ∼ 8 LBGs in the redshift-M_UV plane. For context, we also show recent samples of bright LBGs at similar redshifts from Bouwens et al. (2015, 2016), Calvi et al. (2016), Bowler et al. (2017), Ono et al. (2018), Livermore et al. (2018), and Morishita et al. (2018). Our sample of luminous galaxies is among the most luminous identified at these redshifts, and ∼0.5–1 mag brighter than typical samples selected from CANDELS.

Zoom In Zoom Out Reset image size

In this section we present our measurements of two among the most fundamental observables that the deep near-IR and IRAC observations allow us to investigate, i.e., the spectral slope of the UV-continuum light and the rest-frame u − g color.

The spectral slope of the UV-continuum light is typically parameterized using the so-called UV-continuum slope β (where β is defined such that ${f}_{\lambda }\propto {\lambda }^{\beta }$, Meurer et al. 1999). A common way of deriving the UV-continuum slope is by considering power-law fits to all photometric constraints in the UV continuum (Bouwens et al. 2012; Castellano et al. 2012). Here we take a slightly different approach. First we derive βs for a grid of redshifted Bruzual & Charlot (2003, hereafter BC03) stellar population models with an age of 10 Myr and a range of visual attenuation A_V = 0–2 mag. Then for each individual galaxy we fit the predicted J, H, and K_s-band fluxes to the observations. Uncertainties are derived by randomly scattering the observed fluxes and photometric redshifts by their errors and refitting. This procedure allows us to make full use of the near-IR data and to naturally take into account redshift uncertainties and the Lyman-break entering the J-band at z > 8.5. We caution that, for a small fraction of sources with z > 8.5, βs derived in this way could still be affected by the Lyα emission line shifting into the J-band. We note, however, that observed Lyα equivalent widths of bright z ∼ 8–9 galaxies are modest, 10–30 Å (Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016). As an exercise, we also computed the UV slopes by directly fitting the power law to the flux densities in those bands whose effective wavelength was redder than the redshifted 1300 Å of each object (typically J, H, and K_s). These new estimates (β_phot) resulted in values essentially equal to those from the method we initially applied (median β_phot−β_BC03 ∼ 0.1), although with large scatter for ∼30% of the sources (Δβ ≳ 1). Nonetheless, the large associated uncertainties make the two measurements consistent with each other. However, we believe that the UV slope measurements recovered with the initial method are more robust, because they better model the effects of redshift on the observed flux density of each source.

Figure 10 shows the distribution of UV slopes β and rest-frame u − g colors for the bright z ∼ 8 sample. The z ∼ 8 galaxies span a substantial range in UV spectral slope and color. The large uncertainties, however, suggest that the observed scatter is likely the combination of intrinsic variation and measurement uncertainties. The average slope of the UV continuum, β = −2.2 ± 0.6, is bluer but still consistent with the UV−continuum slopes found for bright −22 < M_UV < −21 galaxies at z = 6 (β = −1.55 ± 0.17) and z ∼ 7 (β = −1.75 ± 0.18) by Bouwens et al. (2014) and suggests a continuing trend toward bluer βs at higher redshifts.

Zoom In Zoom Out Reset image size

Recently, Oesch et al. (2013) analyzed the rest-frame UV and optical properties of a sample of z ∼ 4 LBGs selected from the GOODS-N/S field and Hubble Ultra-Deep Field (HUDF) and spanning a wide range of UV luminosities, M_UV ∼ −18 to ∼ −22 AB. Their J₁₂₅ − [4.5] color (corresponding to approximately rest-frame u − z at z ∼ 4) shows a correlation with the UV slope β (see e.g., their Figure 4), likely driven by dust extinction. The uniform scatter observed at z ∼ 8, then, may suggest rapidly evolving physical mechanisms responsible for the production of dust during the ∼800 Myr between the two epochs.

Recent observational studies have found that the [3.6]–[4.5] color of galaxies depends dramatically on the redshift of the source (Shim et al. 2011; Labbé et al. 2013; Stark et al. 2013; Bowler et al. 2014; Smit et al. 2014, 2015; Faisst et al. 2016; Harikane et al. 2018), with some sources showing extreme colors (Ono et al. 2012; Finkelstein et al. 2013; Laporte et al. 2014, 2015; Faisst et al. 2016; Roberts-Borsani et al. 2016). A number of works have suggested that these extreme colors are likely due to very strong line emission (Labbé et al. 2013; Smit et al. 2014) whereas the intrinsic color of the stellar continua in the absence of emission lines is [3.6]–[4.5] ∼ 0 mag (Labbé et al. 2013; Smit et al. 2014; Rasappu et al. 2016).

At redshift z = 7.0–9.1, the [O iii]+Hβ line emission contributes to the Spitzer/IRAC 4.5 μm band in galaxies, producing red [3.6]–[4.5] colors. Figure 11 shows examples of model colors as a function of redshift for lines with very high equivalent width. Using a small sample of z ∼ 8 galaxies selected from the CANDELS survey, Roberts-Borsani et al. (2016) reported a very red median [3.6]–[4.5] ∼ 0.8 mag color at bright H < 26 mag. Using a simple spectral model, consisting of a flat rest-frame 0.3–0.6 μm continuum in f_ν (i.e., a continuum [3.6]–[4.5] = 0 mag or f_λ ∝ λ⁻²), with the strongest emission lines ([O ii]₃₇₂₇, Hβ, [O iii]_4959,5007, Hα, [N ii]_6548,6583, [S ii]_6716,6730), empirical emission lines ratios from Anders & Alvensleben (2003) for 0.2 Z_⊙ metallicity, they inferred a median [O iii]+Hβ EW of ∼2000 Å. However, the sample of Roberts-Borsani et al. (2016) was very small, and possibly biased because it was compiled from IRAC-selected [3.6]–[4.5] > 0.5 galaxies and galaxies with confirmed Lyα emission. So it is unclear if those results were representative of the general bright z ∼ 8 population.

Zoom In Zoom Out Reset image size

With the UltraVISTA sample and the deep IRAC observations from SPLASH, SEDS, and SMUVS, we have an opportunity to revisit the analysis of Roberts-Borsani et al. (2016) with a larger sample. In Figure 11, we present the [3.6]–[4.5] color distribution for bright z ∼ 8 galaxies from both our study and that of Roberts-Borsani et al. (2016). The [3.6]–[4.5] color distribution spans a range of more 1 mag, with the UltraVISTA sample showing a median [3.6]–[4.5] = 0.62 mag; this color remains unchanged when also combining it with the CANDELS sample.

Adopting the same model of Roberts-Borsani et al. (2016) (see also Smit et al. 2014) and supposing that the 3.6 μm band receives only a negligible contribution from line emission, a [3.6]–[4.5] color of ∼0.6 mag corresponds to an [O iii]+Hβ EW of ∼1500 Å. Such a result is consistent with Labbé et al. (2013) and Smit et al. (2014, 2015), and with the recent estimates of Stefanon et al. (2019, in preparation) and de Barros et al. (2018) based on samples of z ∼ 8 L < L* LBGs selected over the GOODS-N/S fields, which benefit from among the deepest IRAC 3.6 and 4.5 μm observations of the GREATS program (PI: I. Labbé; Labbé et al. 2018, in preparation).

Under the assumption that the extreme IRAC colors are due to nebular emission, our results combined with those from the literature indicate that strong emission lines might be ubiquitous at these redshifts in galaxies spanning ∼3 mag range in luminosity. Nevertheless, significant systematic uncertainties remain depending on the assumed continuum shape and line flux ratios. For example, including the full line list of Anders & Alvensleben (2003) and contribution from the higher-order Balmer lines, assuming a more realistic spectral continuum (e.g., BC03 and scaling emission lines by the flux in hydrogen ionizing photons N_LyC), and allowing for Calzetti et al. (2000) dust, produces a different [3.6]–[4.5] color versus redshift relation by up to 0.2–0.4 mag. Also, emission line ratios, in particular [O iii]₅₀₀₇, depend strongly on metallicity (e.g., Inoue 2011). Considering this, we estimate that simple approximations are probably uncertain by factors of 2–3.

In this section we present our estimates of stellar population parameters for the bright z ∼ 8 galaxies. Measurements were performed with the FAST code (Kriek et al. 2009), adopting Bruzual & Charlot (2003) models for subsolar 0.2Z_⊙ metallicity, a Chabrier (2003) IMF, constant star formation, and the Calzetti et al. (2000) dust law. As discussed above, gaseous emission lines contribute significantly to the integrated broadband fluxes. Given standard BC03 models do not include nebular emission, line and continuum nebular emission were added following the procedure of Salmon et al. (2015) and assuming line flux ratios relative to Hβ from the models calculated by Inoue (2011). The luminosity in Hβ is taken to be proportional to the luminosity in hydrogen ionizing photons N_LyC, assuming ionization-recombination equilibrium (case B). The emission line ratios of Inoue (2011) agree well with the empirical compilations of Anders & Alvensleben (2003), with observations of the local galaxy I Zw 18 (Izotov et al. 1999), and the z = 2.3 galaxy from Erb et al. (2010), in particular for the strongest metal line [O iii]₅₀₀₇. In Table 4 we present the results of our stellar population modeling, specifically the stellar mass, star formation rate, specific star formation rate, age, and extinction, together with the UV₁₆₀₀ absolute magnitude, the UV-continuum slope β, and the rest-frame u − g color for each individual candidate bright z ∼ 8 LBG. A summary of the physical properties is presented in Table 5.

Table 4.  Main Physical Parameters for the Sample of Candidate z ∼ 8 LBGs

ID	MUV	UV Slope β	u − g	$\mathrm{log}({M}_{\star })$	$\mathrm{log}(\mathrm{SFR})$	$\mathrm{log}(\mathrm{sSFR})$	$\mathrm{log}(\mathrm{age})$	AV
	(mag)		(mag)	(M⊙)	(M⊙ yr−1)	(yr−1)	(yr)	(mag)
UVISTA-Y1	−22.48 ± 0.15	$-{1.5}_{-0.7}^{+0.4}$	0.45 ± 0.18	${10.0}_{-0.4}^{+0.9}$	${1.59}_{-9.55}^{+1.02}$	$-{8.4}_{-9.8}^{+1.8}$	${7.30}_{-0.61}^{+1.42}$	${0.9}_{-0.9}^{+0.0}$
UVISTA-Y2	−22.37 ± 0.20	$-{2.6}_{-0.5}^{+0.5}$	0.66 ± 0.21	${9.0}_{-1.2}^{+0.3}$	${1.98}_{-7.57}^{+0.65}$	$-{7.0}_{-8.6}^{+0.8}$	${7.00}_{-0.50}^{+1.80}$	${0.4}_{-0.4}^{+0.3}$
UVISTA-Y3aa	−21.77 ± 0.32	$-{1.5}_{-0.9}^{+0.7}$	0.57 ± 0.39	${9.8}_{-0.3}^{+1.3}$	$-{1.34}_{-7.01}^{+4.06}$	$-{11.1}_{-7.0}^{+4.9}$	${8.00}_{-1.50}^{+0.80}$	${0.0}_{-0.0}^{+1.1}$
UVISTA-Y3ba	−21.23 ± 0.54	$-{3.2}_{-0.0}^{+2.2}$	−0.22 ± 1.82	${8.7}_{-0.0}^{+0.1}$	$-{0.28}_{-0.07}^{+0.00}$	$-{9.0}_{-0.0}^{+0.0}$	${7.40}_{-0.02}^{+0.00}$	${0.0}_{-0.0}^{+0.0}$
UVISTA-Y3ca	−21.37 ± 0.53	$-{2.0}_{-0.0}^{+1.7}$	0.81 ± 0.57	${10.2}_{-0.5}^{+1.6}$	${1.69}_{-28.23}^{+1.86}$	$-{8.5}_{-28.1}^{+2.3}$	${8.70}_{-2.20}^{+0.10}$	${0.8}_{-0.8}^{+1.1}$
UVISTA-Y4	−22.11 ± 0.24	$-{2.7}_{-0.4}^{+0.7}$	0.47 ± 0.26	${9.9}_{-0.2}^{+0.5}$	${1.23}_{-6.97}^{+0.61}$	$-{8.6}_{-7.0}^{+0.9}$	${8.50}_{-1.31}^{+0.30}$	${0.0}_{-0.0}^{+0.7}$
UVISTA-Y5	−22.34 ± 0.24	$-{1.7}_{-0.9}^{+0.8}$	0.63 ± 0.27	${9.0}_{-1.1}^{+0.4}$	${1.99}_{-7.68}^{+0.55}$	$-{7.0}_{-8.6}^{+0.8}$	${7.30}_{-0.80}^{+1.40}$	${0.4}_{-0.4}^{+0.2}$
UVISTA-Y6	−21.92 ± 0.26	$-{1.7}_{-0.8}^{+0.7}$	0.49 ± 0.32	${9.7}_{-0.5}^{+1.1}$	${1.36}_{-12.70}^{+1.33}$	$-{8.4}_{-12.9}^{+2.2}$	${7.30}_{-0.80}^{+1.50}$	${0.9}_{-0.9}^{+0.3}$
UVISTA-Y7	−21.74 ± 0.36	$-{2.0}_{-0.5}^{+0.7}$	⋯b	⋯b	⋯b	⋯b	⋯b	⋯b
UVISTA-Y8	−21.76 ± 0.35	$-{2.8}_{-0.4}^{+0.9}$	0.91 ± 0.45	${8.3}_{-1.4}^{+0.1}$	${1.90}_{-1.41}^{+0.35}$	$-{6.4}_{-2.6}^{+0.2}$	${6.50}_{-0.00}^{+2.29}$	${0.0}_{-0.0}^{+0.5}$
UVISTA-Y9	−21.66 ± 0.34	$-{2.6}_{-0.6}^{+0.9}$	⋯b	⋯b	⋯b	⋯b	⋯b	⋯b
UVISTA-Y10	−21.89 ± 0.31	$-{2.2}_{-0.7}^{+1.1}$	0.56 ± 0.39	${8.3}_{-1.4}^{+0.0}$	${1.80}_{-4.34}^{+0.47}$	$-{6.5}_{-5.6}^{+0.3}$	${6.70}_{-0.20}^{+2.10}$	${0.0}_{-0.0}^{+0.5}$
UVISTA-Y11	−22.04 ± 0.26	$-{1.8}_{-1.3}^{+0.5}$	0.67 ± 0.30	${8.7}_{-1.2}^{+0.4}$	${1.76}_{-7.62}^{+0.60}$	$-{7.0}_{-8.6}^{+0.8}$	${7.30}_{-0.80}^{+1.44}$	${0.3}_{-0.3}^{+0.2}$
UVISTA-Y12	−21.66 ± 0.40	$-{2.1}_{-0.8}^{+1.3}$	0.22 ± 0.64	${9.1}_{-0.4}^{+0.9}$	${0.17}_{-2.88}^{+2.22}$	$-{9.0}_{-3.1}^{+2.8}$	${7.40}_{-0.90}^{+1.30}$	${0.2}_{-0.2}^{+0.3}$
UVISTA-Y13	−21.39 ± 0.42	$-{1.0}_{-0.7}^{+0.7}$	0.50 ± 0.51	${9.8}_{-0.3}^{+1.3}$	${0.70}_{-9.08}^{+1.82}$	$-{9.1}_{-9.1}^{+2.8}$	${7.50}_{-0.96}^{+1.28}$	${0.8}_{-0.8}^{+0.3}$
UVISTA-Y14	−21.44 ± 0.40	$-{3.0}_{-0.1}^{+1.7}$	0.63 ± 0.54	${9.3}_{-0.4}^{+1.2}$	${0.52}_{-9.14}^{+2.03}$	$-{8.8}_{-9.3}^{+2.6}$	${8.20}_{-1.70}^{+0.63}$	${0.0}_{-0.0}^{+0.8}$
UVISTA-Y15	−21.50 ± 0.37	$-{2.7}_{-0.4}^{+0.9}$	−0.04 ± 0.68	${8.8}_{-0.0}^{+0.2}$	$-{0.16}_{-0.02}^{+0.41}$	$-{9.0}_{-0.0}^{+0.6}$	${7.40}_{-0.10}^{+0.01}$	${0.0}_{-0.0}^{+0.0}$
UVISTA-Y16	−21.80 ± 0.33	$-{2.2}_{-0.8}^{+0.8}$	0.39 ± 0.40	${8.6}_{-0.9}^{+0.3}$	${1.62}_{-0.98}^{+0.56}$	$-{7.0}_{-1.8}^{+0.8}$	${7.30}_{-0.80}^{+1.50}$	${0.1}_{-0.1}^{+0.4}$

Notes.

^aThese three candidate LBGs were originally identified as a single source, successively deblended using data from the COSMOS/DASH program (see Section 5.1 and Figure 8). When we do not deblend the source, we obtain M_UV = −22.00 ± 0.16 mag, β = −1.8 ± 0.7, u − g = 0.58 ± 0.16 mag, $\mathrm{log}({M}_{\star }/{M}_{\odot })={9.9}_{-0.3}^{+0.6}$, $\mathrm{log}(\mathrm{SFR}/{M}_{\odot }/{\mathrm{yr}}^{-1})={1.63}_{-3.77}^{+0.38}$, $\mathrm{log}(\mathrm{sSFR}/{\mathrm{yr}}^{-1})=-{8.2}_{-3.8}^{+0.9}$, $\mathrm{log}(\mathrm{age}/\mathrm{yr})={8.20}_{-1.16}^{+0.60}$, and ${A}_{V}={0.5}_{-0.5}^{+0.5}$ mag. ^bAfter visual inspection, the neighbor-cleaned image stamps in the IRAC 3.6 and 4.5 μm bands showed non-negligible residuals that likely systematically affected our estimates. Photometric redshifts resulted to be robust against the exclusion of the flux densities in these two bands, but stellar population parameters heavily rely on the IRAC colors. Because of the unreliability of the IRAC flux density estimates for these objects, we discard their physical parameters.

Download table as:  ASCIITypeset image

Table 5.  Observed and Rest-frame Properties for Candidate z ∼ 8 Galaxies Identified in the UltraVISTA DR3 Observations

Quantity	25%	Median	75%	25% Uncertainties	Median Uncertainties	75% Uncertainties
zphot	8.05	8.40	8.62	+0.60/−0.61	+0.69/−0.73	+0.96/−1.15
MUV [mag]	−22.0	−21.8	−21.6	±0.3	±0.3	±0.4
UV β	−2.68	−2.17	−1.73	+0.70/−0.40	+0.77/−0.65	+1.01/−0.79
(u−g)rest [mag]	0.42	0.53	0.65	±0.29	±0.39	±0.55
$\mathrm{log}({M}_{\star }/{M}_{\odot })$	8.71	9.07	9.76	+0.32/−0.24	+0.46/−0.44	+1.14/−1.15
M⋆/LUV [M⊙/L⊙]	0.005	0.010	0.044	+0.008/−0.004	+0.015/−0.014	+0.034/−0.098
M⋆/Lu [M⊙/L⊙]	0.013	0.048	0.101	+0.026/−0.010	+0.038/−0.056	+0.098/−0.244
M⋆/Lg [M⊙/L⊙]	0.017	0.064	0.133	+0.026/−0.007	+0.043/−0.089	+0.078/−0.188
log(SFR/M⋆/yr−1)	0.3	1.5	1.8	+0.5/−2.1	+0.6/−7.3	+1.8/−9.1
$\mathrm{log}(\mathrm{sSFR}/{\mathrm{yr}}^{-1})$	−9.0	−8.4	−7.0	+0.7/−2.9	+0.9/−7.8	+2.5/−9.2
$\mathrm{log}(\mathrm{age}/\mathrm{yr})$	7.30	7.35	7.75	+0.47/−0.35	+1.35/−0.80	+1.50/−1.13
AV [mag]	0.00	0.15	0.60	+0.22/−0.00	+0.32/−0.15	+0.56/−0.60

Note. Estimates of z_phot, M_UV, and L_X were obtained from EAzY (see Section 4); M_⋆, SFR, sSFR, age, and A_V were measured with FAST (see Section 6.4); and the UV-continuum slope β were measured following the procedure described in Section 6.2. The last two columns present the first and third quartiles of uncertainties, respectively.

Download table as:  ASCIITypeset image

As we already introduced in Section 4, the neighbor-cleaned IRAC 3.6 and 4.5 μm band image sections for two sources (UVISTA-Y7 and UVISTA-Y9) presented residuals that might be systematically affecting our estimates of stellar population parameters (see Figure 6). We therefore recomputed the redshift likelihood distributions for these two sources after excluding the IRAC flux densities. The photometric redshifts we derived were consistent with the estimates obtained adopting the full set of measurements. However, the stellar population parameters heavily rely on the IRAC colors because at z ∼ 8 these probe the rest-frame optical redward of the Balmer break and the emission line properties, both affecting their age and the stellar mass measurements. As a result, the physical parameters for the two sources have not been included in Table 4 or figures presenting these parameters (i.e., Figures 10–13).

Zoom In Zoom Out Reset image size

In Section 6.3 we showed that our sample is characterized by extreme [3.6]–[4.5] ∼ 0.6 mag colors, likely the result of strong [O iii]+Hβ emission entering the 4.5 μm band. A number of studies have shown that nebular emission can systematically bias stellar mass estimates (e.g., Stark et al. 2013). Figure 12 compares the best-fit stellar masses to those derived with the standard BC03 models without emission lines for our sample. Those masses are higher by ∼0.4 dex on average (scatter ∼0.6 dex), with individual galaxies differing by up to 1 dex. This is consistent with Labbé et al. (2013), who estimate that z ∼ 7–8 galaxies' average stellar masses decrease by ∼0.5 dex if the contributions of emission lines to their broadband fluxes are accounted for. However, the discrepancy appears to be related not only to the strong contribution of [O iii]₅₀₀₇ to the 4.5 μm band. Indeed, if we refit the galaxies with the standard BC03 models (without emission lines) while omitting the flux in the 4.5 μm band, the offset is marginally reduced to 0.23 dex (scatter 0.43 dex) compared to the BC03 and emission lines fit to all bands. This residual offset is likely due to the effect of nebular emission (mainly [O ii]₃₇₂₇) characteristic of young stellar populations, which still substantially contaminates the 3.6 μm band. This result stresses once more the importance of accounting for nebular emission in estimating the physical parameters of z ≳ 8 galaxies.

The typical estimated stellar masses for bright sources in our z ∼ 8 selection (see Table 5) are ${10}^{{9.1}_{-0.4}^{+0.5}}$ M_⊙, with the SFRs of ${32}_{-32}^{+44}{M}_{\odot }$ yr⁻¹, specific SFR of ${4}_{-4}^{+8}$ Gyr⁻¹, stellar ages of $\sim {22}_{-22}^{+69}$ Myr, and low dust content ${A}_{V}={0.15}_{-0.15}^{+0.30}$ mag. As evident from Table 5, individual galaxies shows a broad range in each of these properties, with interquartile masses, ages, and specific star formation rates spanning ∼1 dex.

In Figure 13 we compare the rest-frame properties with the best-fit stellar mass-to-light ratios for luminosities in the rest-frame UV₁₆₀₀ and rest-frame g-band. These quantities are not completely independent, because both are derived from the same photometry, but provide useful insights in how color relates to stellar mass. Overall, the mass-to-light ratios are very low, as expected for very young stellar ages (<100 Myr), but span quite a wide range, between 0.1 and 0.01M_⊙/L_⊙.

Zoom In Zoom Out Reset image size

We find a positive although marginal correlation of the M_⋆/L_UV,1600 with the UV slope for our z ∼ 8 sample, because it could be expected from older and/or dustier stellar populations characterized by redder UV slope (e.g., Bouwens et al. 2014).

A number of works have shown that at low redshift there exists a tight relation between rest-frame optical colors and M_⋆/L ratios, such that redder galaxies exhibit higher M_⋆/L, and that this empirical relation is not sensitive to details of the stellar population modeling (e.g., Bell & de Jong 2001). This relation appears to hold even at intermediate redshifts z ∼ 2 (e.g., Szomoru et al. 2013). Remarkably, in contrast to the situation at low-redshift, redder rest-frame u − g colors of the z ∼ 8 sample do not correspond to higher M_⋆/L. Instead, the optically reddest galaxies tend to have the lowest M_⋆/L. This likely reflects the effect of strong emission lines in the g-band. The fact that age and dust have very different effects on the colors of the high-redshift galaxies studied here probably also explains the lack of correlation between β and u − g in Figure 10.

In this section we present our measurements of the UV LF based on the sample presented in this work. Our main result is the UV LF at z ∼ 8 based on the sample of Y-band dropouts (i.e., considering UVISTA-Y3 as three independent sources) presented in Section 5.2. However, because some objects have a nominal photometric redshift z_phot ∼ 9, we also explored the contribution to the UV LF at z ∼ 9 from the five sources with z_phot ≥ 8.6 (namely UVISTA-Y3a, UVISTA-Y3b, UVISTA-Y3c, UVISTA-Y11, and UVISTA-Y12). Because the nominal photometric redshift of UVISTA-Y5 is z_phot = 8.596, this object was initially excluded by our redshift selection criterion. Considering the very marginal difference of its photo-z from the selection threshold, we also forced its inclusion into the sample adopted for the estimate of the z ∼ 9 LF, bringing to six the total number of sources used for the z ∼ 9 LF.

To infer the volume densities of the galaxies, we first estimate the detection completeness and selection function through simulations. Following Bouwens et al. (2015), we generated catalogs of mock sources with realistic sizes and morphologies by randomly selecting images of z ∼ 4 galaxies from the Hubble Ultra Deep Field (Beckwith et al. 2006; Illingworth et al. 2013) as templates. The images were scaled to account for the change in angular diameter distance with redshift and for evolution of galaxy sizes at fixed luminosity ∝(1 + z)⁻¹ (e.g., Oesch et al. 2010; Ono et al. 2013; Holwerda et al. 2015; Shibuya et al. 2015). The template images are then inserted into the observed images, assigning colors expected for star-forming galaxies in the range 6 < z < 11. The colors were based on a UV-continuum slope distribution of β = −1.8 ± 0.3 to match the measurements for luminous 6 < z < 8 galaxies (Bouwens et al. 2012, 2014; Finkelstein et al. 2012; Rogers et al. 2014). The simulations include the full suite of HST, ground-based, and Spitzer/IRAC images. For the ground-based and Spitzer/IRAC data the mock sources were convolved with appropriate kernels to match the lower-resolution PSF. To simulate IRAC colors we assume a continuum flat in f_ν and strong emission lines with fixed rest-frame EW(Hα+[N ii]+[S ii]) = 300 Å and rest-frame EW([O iii]+Hβ) = 500 Å consistent with the results of Labbé et al. (2013), Stark et al. (2013), Smit et al. (2014, 2015), and Rasappu et al. (2016). We included the effect of other nebular lines following the recipe of Anders & Alvensleben (2003) for subsolar metallicity.

The same detection and selection criteria as described in Section 4 were then applied to the simulated images to calculate the completeness as a function of recovered magnitude and the selection as a function of magnitude and redshift (see Figure 8 of Stefanon et al. 2017b for the selection functions over the UltraVISTA deep and ultradeep stripes).

The total selection volume over our UltraVISTA area for galaxies with H ∼ 24.0–24.5 mag and 24.5–25.0 is 5.3 × 10⁶ Mpc³ and 2.6 × 10⁶ Mpc³, respectively.

We estimate constraints on the bright end of the UV LF, adopting the V_max formalism of Avni & Bahcall (1980) in 0.5 mag bins, optimizing the range in UV luminosities of the sample. Following Moster et al. (2011), we increase by 24% the Poisson uncertainties to account for cosmic variance. The resulting z ∼ 8 LF is shown in the top panel of Figure 14, and the corresponding number densities are listed in Table 6. In our measuring, we only included sources more luminous than M_UV ≤ −21.3 mag, for a total of 17 sources, excluding UVISTA-Y3b due to its extremely low luminosity, which makes the estimate of the completeness at that luminosity uncertain. Nonetheless, we stress that the volume density we derive in the faintest luminosity bin is likely a lower limit, as the actual incompleteness may be larger than what we estimate. Considering that six sources in our sample are characterized by redshifts z ≥ 8.6 (after forcing the inclusion of UVISTA-Y5), we considered these galaxies to belong to the z ∼ 9 redshift bin and computed the associated number densities accordingly. The resulting z ∼ 9 LF is presented in the bottom panel of Figure 14 and in Table 6.

Zoom In Zoom Out Reset image size

Table 6.  V_max Determinations of the UV LF

MUV	ϕ
(mag)	(×10−6 mag−1 Mpc−3)
z ∼ 8
−22.55	${0.76}_{-0.41}^{+0.74}$
−22.05	${1.38}_{-0.66}^{+1.09}$
−21.55a	${4.87}_{-1.41}^{+2.01}$
z ∼ 9
−22.35b	${0.43}_{-0.36}^{+0.99}$
−22.00	${0.43}_{-0.36}^{+0.98}$
−21.60	${1.14}_{-0.73}^{+1.50}$
−21.20a	${1.64}_{-1.06}^{+2.16}$

Notes.

^aThis luminosity bin includes sources from the deblending of UVISTA-Y3, which fall below our nominal detection threshold. The sample in this luminosity bin is therefore likely incomplete. ^bThe volume density in this luminosity bin was obtained forcing UVISTA-Y5 into the sample of galaxies at z_phot ≥ 8.6 (i.e., our z ∼ 9 LBG sample). Its nominal z_phot = 8.596 would exclude it from the sample of z ∼ 9 sources when the redshift selection criteria is strictly enforced; however, considering the very small difference with the z_phot threshold, here we include it for completeness.

Download table as:  ASCIITypeset image

In Figure 14 we also compare our LF estimates with other recent estimates of the bright end of the LF from empty field searches at z ∼ 8 (Bradley et al. 2012; McLure et al. 2013; Schenker et al. 2013; Schmidt et al. 2014; Bouwens et al. 2015; Finkelstein et al. 2015; Roberts-Borsani et al. 2016; Stefanon et al. 2017b; Bridge et al. 2019) and z ∼ 9 (Oesch et al. 2013; Bouwens et al. 2016; Calvi et al. 2016; McLeod et al. 2016; Ishigaki et al. 2018; Livermore et al. 2018; Morishita et al. 2018). The volume density of z ∼ 8 LBGs probed here corresponds to a luminosity range that exhibits only a modest overlap with earlier LF studies (i.e., Bouwens et al. 2010, 2011; McLure et al. 2013; Schenker et al. 2013; Schmidt et al. 2014; Finkelstein et al. 2015), where essentially all z ∼ 8 candidates have apparent magnitudes fainter than H ∼ 25.5 mag. Nonetheless, our luminosity regime overlaps with the widest-area searches available to date from the CANDELS fields (Bouwens et al. 2015 and Roberts-Borsani et al. 2016, which includes the spectroscopically confirmed z ∼ 8 LBGs of Oesch et al. 2015b and Zitrin et al. 2015) and from the BoRG program (Trenti et al. 2011; Calvi et al. 2016; Livermore et al. 2018; Morishita et al. 2018; Bridge et al. 2019).

Perhaps quite unsurprisingly, the new estimate of the z ∼ 8 LF is consistent with the previous measurement at M_UV ≲ −22 mag of Stefanon et al. (2017b), based on a partly different analysis of the six among the brightest sources presented in this work (UVISTA-Y1 through UVISTA-Y6), and where we also considered HST/WFC3 imaging for three of them from one of our HST programs. The availability of HST/WFC3 DASH data allowed us to ascertain that UVISTA-Y3 is likely a triple system of fainter (∼L*) LBGs. However, the revised analysis performed for the current work showed that one of the sources previously considered to be at z_phot ∼ 7.5 is actually at z_phot ≳ 8, thus increasing its luminosity and balancing the final volume density.

For M_UV ≲ −22 mag sources, our new results are also consistent with the upper limits of Bradley et al. (2012), Bouwens et al. (2015), Finkelstein et al. (2015), and Roberts-Borsani et al. (2016); our measurements are in excess of what is expected, extrapolating the Bouwens et al. (2015) results to brighter magnitudes by a factor of 8, but are nevertheless consistent within 2σ.

At M_UV ≳ −22 mag our new estimates are consistent with the volume densities of bright LBGs over the CANDELS fields reported by McLure et al. (2013), Bouwens et al. (2015), Finkelstein et al. (2015), and Roberts-Borsani et al. (2016), and with the measurements of Bradley et al. (2012), Schenker et al. (2013), and Schmidt et al. (2014) from the BoRG program (Trenti et al. 2011; Yan et al. 2011). Recently, Bridge et al. (2019) presented the z ∼ 8 LF from eight M_UV ≳ − 22 mag sources identified over BoRG fields for which Spitzer/IRAC data were collected in the 3.6 and 4.5 μm bands. The associated volume density is ∼5× higher than what we estimate for our sample, and their measurements are only consistent at ∼3σ. However, the steepness of the LF at the bright end significantly increases the challenges in comparing volume density estimates due to the sensitive dependence on the precise luminosity range probed in different studies. Furthermore, this discrepancy could in part be explained by the different median cosmic times probed by the two samples, considering that the median redshift of the Bridge et al. (2019) sample, z_phot,med = 7.76, is lower than the median redshift of our sample (z_phot ∼ 8.4).

In the lower panel of Figure 14 we present our estimates of the z ∼ 9 LF. Here we mark with an open symbol the point corresponding to the faintest bin of luminosity, because our selection in that luminosity range is likely very incomplete.

At M_UV ≲ −22 mag our new bright z ∼ 9 results are consistent with the upper limits of Bouwens et al. (2016) from CANDELS and of Ishigaki et al. (2018) from the HFF initiative (Lotz et al. 2017). Our measurement at M_UV ∼ −22 mag is consistent at 1σ with the measurement of Morishita et al. (2018) based on BoRG observations partly supported by Spitzer/IRAC observations, while it is consistent with that of Calvi et al. (2016) at ∼2σ, our density being ∼9× lower than the corresponding measurement of Calvi et al. (2016). One possible explanation for differences between our results and those of Calvi et al. (2016) would be if the Calvi et al. (2016) samples suffer from significant contamination. This is especially a concern because few of candidate z ∼ 9 sources have available Spitzer/IRAC or deep Y₀₉₈ imaging to aid in source selection. In fact, Livermore et al. (2018) find that one especially bright candidate reported by Calvi et al. (2016) appeared to be clearly a low-redshift candidate after further examination.

In the same panel we also plot a double power law (DPL) that we evolved to z ∼ 9, applying the relations of Bouwens et al. (2016) to the DPL found at z ∼ 7 by Bowler et al. (2017, see also Stefanon et al. 2017b). Indeed, the excess in number density we observe for M_UV < −22 mag at z ∼ 9 (introduced by forcing UVISTA-Y5 into the z ∼ 9 sample) seems to be better described by the double power law. We remark, however, that the still large uncertainties do not allow us to fully remove the degeneracy on the shape of the LF at z ∼ 9, which instead needs larger samples. We will discuss the shape of the z ∼ 8 LF in more detail in Section 6.7.

The bright candidates found over UltraVISTA alone are not sufficient to constrain the overall shape of the UV LF due to lack of dynamic range. In the case of a Schechter (1976) function where the shape is determined by the faint-end slope α and turn over magnitude M*, both bright and faint objects are needed to constrain α and M*. The similar redshift distributions expected for bright galaxies selected by our z ∼ 8 criteria and those selected in the fainter Bouwens et al. (2015) samples (see Figure 2) make it possible to combine our z ∼ 8 LF with the corresponding estimates of Bouwens et al. (2015), based on the CANDELS, HUDF09, HUDF12, ERS, and BoRG/HIPPIES programs.

The combined step-wise determination of the UV LF at z ∼ 8 is presented in Figure 15. We determined the Schechter function parameters M*, α, and ϕ* minimizing the χ², and obtaining $\mathrm{log}({\phi }^{* })=-{3.99}_{-0.37}^{+0.29}$, ${M}_{1600}^{* }=-{20.95}_{-0.35}^{+0.30}$ mag, and $\alpha =-{2.15}_{-0.19}^{+0.20}$. The 68% and 95% confidence level contours are presented in Figure 16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our sample of bright z ∼ 8 LBGs make the characteristic luminosity M* brighter by ∼0.5 mag compared to the most recent estimates of Bouwens et al. (2015), even though this result is significant only at ∼1.2σ, while the faint-end slope α is consistent at 1σ.

In Figure 16 we also compare our estimated Schechter parameters to their evolution over a wide range of redshift, 4 ≲ z ≲ 10 from Bouwens et al. (2015). Our result confirms the picture of marginal evolution of M* for z ≳ 3–4, but significant evolution of α and ϕ*. This conclusion was first drawn by Bouwens et al. (2015) using LF results from z ∼ 7 to z ∼ 4 (see also Finkelstein et al. 2015), although they are in modest tension with the results of Bowler et al. (2015), who suggest an evolution of dM*/dz ∼ 0.2 from z ∼ 7 to z ∼ 5.

One significant area of exploration over the last few years has regarded the form on the UV LF at the bright end. In particular, there has been interest in determining whether the UV LF shows more of an exponential cut-off at the bright end or a power-law-like cut-off. The higher number densities implied by a power-law-like form might indicate that the impact of either feedback or dust is less important at high redshifts than it is at later cosmic times. Successfully distinguishing a power-law-like form for the bright end of the LF from a sharper exponential-like cut-off is challenging, because it requires very tight constraints on the bright end of the LF and hence substantial volumes for progress.

The simplest functional form to use in fitting the UV LF is a power law, which can be useful when very wide-area constraints are not available for fitting the bright end. One of the earliest considerations of a power-law form in fitting the UV LF at z > 6 was by Bouwens et al. (2011), and it was shown that such a functional form satisfactorily fit all constraints on the z ∼ 8 LF from HST available at the time (Figure 9 from that work).

Here we consider three functional forms that can potentially be adopted to describe the number density of galaxies at z ∼ 8: a single power law, a DPL, and the Schechter (1976) form. The parameterization for a double power law is as follows (see also Bowler et al. 2012; Ono et al. 2018):

$$\begin{eqnarray*}&&\phi (M)=\displaystyle \frac{{\phi }^{* }}{{10}^{0.4(\alpha +1)(M-{M}^{* })}+{10}^{0.4(\beta +1)(M-{M}^{* })}}\end{eqnarray*}$$

where α and β are the faint-end and bright-end slopes, respectively, M* is the transition luminosity between the two power-law regimes, and ϕ* is the normalization.

A quick inspection of Figure 15 suggests already that the UV LF at z ∼ 8 cannot be well represented by a power-law form. Indeed, a χ² test, as previously adopted by e.g., Bowler et al. (2012, 2015, 2017) at z ∼ 7, results in reduced χ², ${\chi }_{\nu }^{2}\,=3.5,1.04$, and 1.05 for the single power law, DPL, and Schechter functional form, respectively. The double-power-law parameters are α = −1.92 ± 0.50, β = −3.78 ± 0.48, M* = −20.04 ± 1.00 mag, and ${\phi }^{* }={3.88}_{-3.88}^{+5.80}\times {10}^{-4}$ Mpc⁻³ mag⁻¹.

The above results suggest that we cannot yet properly distinguish between a Schechter and a double-power-law form, a result that might be driven by the higher volume density we measured in the brightest absolute magnitude bin. Nevertheless, this result is in line with recent UV LF estimates at z ≲ 7 from large area surveys (UltraVISTA DR2—Bowler et al. 2014, 2015, HSC Survey—Ono et al. 2018), which found an excess in the volume densities of z ∼ 4–7 galaxies for L > L* compared to the Schechter exponential decline.

Even though our favored interpretation consists in considering UVISTA-Y3 composed by three independent sources, in Appendix C, for completeness, we also present the V_max estimates obtained from the blended UVISTA-Y3. We note, however, that these new estimates do not change significantly from those presented in this section.

In Section 5.5 we identified four sources whose flux was likely amplified by massive foreground galaxies. Indeed, recent studies have found that gravitational lensing magnification could explain, at least in part, the excess in number density observed at the bright end of z ∼ 4−7 UV LF (see e.g., Ono et al. 2018). Correcting the apparent magnitude of the four impacted sources for the estimated magnitude and re-deriving the Schechter parameters from the z ∼ 8 constraints, we find M* = −20.81 mag, ∼0.14 mag fainter than with no correction (see Appendix D for details), possibly indicating that lensing is likely playing a modest role in shaping the bright end (M_UV ≲ −22.5 mag) of the z ∼ 8 LF. Because the impact of lensing amplification is uncertain and also model dependent, we follow Bowler et al. (2015) in ignoring the impact for our fiducial determinations.

An alternative way of making sense of the overall shape of the UV LF is to compare it to the HMF. To this aim, we scaled the HMF by a fixed M_halo/L_UV ratio to match the UV LF. We present the result in Figure 15, adopting the Sheth et al. (2001) HMF generated by HMFcalc¹⁶ (Murray et al. 2013), assuming a Planck Collaboration et al. (2016) cosmology, with Ω_m = 0.2678, Ω_b = 0.049, H₀ = 67.04, n_s = 0.962, and σ₈ = 0.8347.

The scaled HMF looks similar to the UV LF at z ∼ 8. We observe only a small Δα ∼ 0.4 difference in the faint-end slope. We also observe a slight difference in the effective slope at the bright end when a Schechter form is considered (Δβ ∼ 0.2–0.3). Interestingly, the bright end of the double power law overlaps with the HMF for M_UV ≲ −22 mag, and might suggest different feedback efficiencies (see also Bowler et al. 2014; Ono et al. 2018). However, as we concluded earlier in this section, our results do not allow us to ascertain whether the cut-off we observe is exponential in form or has a more power-law-like form.

We will not conduct a similar quantitative assessment of the shape of the UV LF at z ∼ 9, due to the challenges in determining the total number of bright z ∼ 9 galaxies over UltraVISTA. Clearly, if any significant number of the candidates do prove to be bona fide z ∼ 9 galaxies, they would favor more of a power-law form to the bright end of the z ∼ 9 LF.

The evolution of the UV luminosity density with cosmic time provides an insight into the rate at which galaxies are building up and how this rate might depend on cosmic time or galaxy/halo mass.

In Figure 17, we present the UV luminosity densities we infer to the effective faint-end limit of our present search (i.e., ∼−21.5 mag) at z ∼ 8 obtained assuming a Schechter form, given a power law at the bright end would imply an infinite luminosity density. We note, however, that we obtain identical results if we compute the luminosity density over M_UV = [−23, −21.5] mag with both the Schechter and DPL parameterization. The luminosity density constraint we find at z ∼ 8 is ${10}^{{23.87}_{-0.68}^{+0.58}}$ ergs/s/Hz/Mpc³ (68% c.l.).

Zoom In Zoom Out Reset image size

In Figure 17 we also present the luminosity density determinations and uncertainties that we derive integrating to M_UV = −21.5 mag the recent LF results of Bouwens et al. (2015) at z ∼ 4, z ∼ 5, z ∼ 6, z ∼ 7, and z ∼ 8; the results of Bouwens et al. (2016) at z ∼ 9; and the results at z ∼ 10 from Oesch et al. (2018), together with our best-fit constraints and 1σ uncertainties on the evolution of the luminosity density with redshift, assuming that the logarithm of the luminosity density decreases linearly with increasing the redshift.

Our estimate is consistent with that derived from the z ∼ 8 LF of Bouwens et al. (2015) and with the upper limit at z ∼ 9 based on Bouwens et al. (2016) LF. However, the luminosity densities recovered from the LF of Oesch et al. (2018) and Bouwens et al. (2015) at z ∼ 10 are ∼1.7σ and ∼1.1σ respectively, lower than the extrapolation of the linear relation discussed above, hinting at a potentially even more rapid buildup of the brightest galaxies in the first ∼500 Myr of cosmic history (see e.g., the extensive discussion in Oesch et al. 2018).