<figure><img src="image_1.png" /></figure>

# The Brightest $z \gtrsim 8$ Galaxies over the COSMOS UltraVISTA Field

Mauro Stefanon¹, Ivo Labbé², Rychard J. Bouwens¹, Pascal Oesch³, Matthew L. N. Ashby⁴, Karina I. Caputi⁵,⁶, and Henry J. McCracken¹¹

¹Leiden Observatory, Leiden University, NL-2300 RA Leiden, The Netherlands; stefanon@strw.leidenuniv.nl

²Centre for Astrophysics and SuperComputing, Swinburne, University of Technology, Hawthorn, Victoria, 3122, Australia

³Observatoire de Genève, 51 Ch. des Maillettes, 1290 Versoux, Switzerland

⁴Center for Astrophysics—Harvard & Smithsonian 60 Garden St., Cambridge, MA 02138, USA

⁵Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700AV Groningen, The Netherlands

⁶Cosmic Dawn Center (DAWN), Denmark

⁷Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, 2100 Copenhagen Ø, Denmark

⁸UCO/Lick Observatory, University of California, Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA

⁹Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France

¹⁰Department of Physics and Astronomy, Tufts University, Medford, MA 02155, USA

¹¹Sorbonne Université, CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France

¹²York University, 4700 Keele St., Toronto, ON, M3J 1P3, Canada

¹³Astronomy Department, Yale University, 52 Hillhouse Ave, New Haven, CT 06511, USA

Received 2019 February 27; revised 2019 July 19; accepted 2019 July 30; published 2019 September 24

## Abstract

We present 16 new ultrabright $H_{AB} \lesssim 25$ galaxy candidates at $z \sim 8$ identified over the COSMOS/UltraVISTA field. The new search takes advantage of the deepest-available ground-based optical and near-infrared observations, including the DR3 release of UltraVISTA and full-depth *Spitzer*/IRAC observations from the SMUVS and SPLASH programs. Candidates are selected using Lyman-break color criteria, combined with strict optical non-detection and SED-fitting criteria, designed to minimize contamination by low-redshift galaxies and low-mass stars. *HST*/WFC3 coverage from the DASH program reveals that one source evident in our ground-based near-IR data has significant substructure and may actually correspond to 3 separate $z \sim 8$ objects, resulting in a total sample of 18 galaxies, 10 of which seem to be fairly robust (with a >97% probability of being at $z > 7$). The UV-continuum slope $\beta$ for the bright $z \sim 8$ sample is $\beta = -2.2 \pm 0.6$, bluer but still consistent with that of similarly bright galaxies at $z \sim 6$ ($\beta = -1.55 \pm 0.17$) and $z \sim 7$ ($\beta = -1.75 \pm 0.18$). Their typical stellar masses are $10^{9.1^{+0.5}_{-0.4}} M_{\odot}$, with the SFRs of $32^{+44}_{-32} M_{\odot}$ yr$^{-1}$, specific SFR of $4^{+8}_{-4}$ Gyr$^{-1}$, stellar ages of $\sim 22^{+69}_{-22}$ Myr, and low dust content $A_V = 0.15^{+0.30}_{-0.15}$ mag. Using this sample we constrain the bright end of the $z \sim 8$ UV luminosity function. When combined with recent empty field luminosity function estimates at similar redshifts, the resulting $z \sim 8$ luminosity function can be equally well represented by either a Schechter or a double-power-law form. Assuming a Schechter parameterization, the best-fit characteristic magnitude is $M^* = -20.95^{+0.30}_{-0.35}$ mag with a very steep faint-end slope $\alpha = -2.15^{+0.20}_{-0.19}$. These new candidates include some of the brightest objects found at these redshifts, 0.5–1.0 magnitude brighter than those found over CANDELS, and providing excellent targets for spectroscopic and longer-wavelength follow-up studies.

**Key words:** galaxies: evolution – galaxies: formation – galaxies: high-redshift – galaxies: luminosity function, mass function

## 1. Introduction

The confirmation and characterization of galaxy candidates within the cosmic reionization epoch has been a major challenge for observational extragalactic astronomy for the last few years. The exceptional sensitivity offered by the Wide Field Camera 3 Infrared (WFC3/IR) instrument on board the *Hubble Space Telescope* (HST), combined with efficient photometric selection techniques, has enabled the identification of $\gtrsim 700$ faint galaxy candidates at $z = 7–11$ (e.g., Bouwens et al. 2011, 2015; Oesch et al. 2012, 2014, 2016, 2018; McLure et al. 2013; Schenker et al. 2013; Schmidt et al. 2014; Finkelstein et al. 2015). These high-redshift galaxy samples have provided a powerful way to investigate the buildup and evolution of galaxies by imposing new constraints on the evolution of their rest-frame ultraviolet (UV) luminosity functions (LFs) and integrated star formation rate density (SFRD—but see also e.g., Tanvir et al. 2012; McGuire et al. 2016 for a complementary approach using gamma-ray bursts).

The redshift range of $z \sim 8–10$ is of particular interest: a number of works suggest a rapid decline of the SFRD from $z \sim 8$ to $z \sim 10$ (see e.g., Oesch et al. 2012, 2014, 2015a, 2018; Ellis et al. 2013; Bouwens et al. 2015—but see e.g., McLeod et al. 2015, 2016). A key question is therefore whether the faint galaxies emit enough ionizing photons to reionize the universe at $z \gtrsim 7$ (e.g., Bolton & Haehnelt 2007; Oesch et al. 2009; Robertson et al. 2010; Bouwens et al. 2011, 2015; Shull et al. 2012; Finkelstein et al. 2015; Tanvir et al. 2019).

Answering the above question requires estimating the faint-end slope of the UV LF during the reionization epoch. For a Schechter (1976) parameterization of the LF, because of the correlation between the characteristic luminosity and the faint-end slope, constraining the bright end of the LF (e.g., through searches in shallow wide-field surveys) will also improve the estimates at the faint end (e.g., Bouwens et al. 2008). Furthermore, identifying bright Lyman-break galaxies (LBGs)
will help determine whether the LF has an exponential cut-off (with relatively few luminous galaxies, as has been established at $z < 7$) or is featureless like a power law (as suggested by a recent works—e.g., Bowler et al. 2015, 2017; Ono et al. 2018). Finally, measurements of the bright end encode crucial information about early galaxies, including the effects of dust, star formation feedback, and the duty cycle of galaxies. The evolution of the bright end therefore provides strong tests for models of galaxy evolution at these redshifts (e.g., Finlator et al. 2011; Jaacks et al. 2012; Mason et al. 2015; Trac et al. 2015; Mashian et al. 2016; Waters et al. 2016).

Bright $z \gtrsim 8$ candidate LBGs are also important targets for spectroscopic follow-up and in preparation for the James Webb Space Telescope. Spectroscopic confirmation is vital to test the validity of the photometric selection techniques and to identify potential contaminant populations at lower redshift, given the physical conditions at such early times are potentially very different than at present, increasing the uncertainty in photometric redshift determinations. When galaxies are confirmed, spectroscopy enables the study of UV spectral features (e.g., Stark et al. 2015a, 2015b, 2017) and improves estimates of stellar mass and star formation rate. However, spectroscopic confirmation has been very challenging so far, with fewer than expected (e.g., Stark et al. 2011) normal galaxies with robust redshift measurements at $z > 7$ (e.g., Pentericci et al. 2011; Vanzella et al. 2011; Ono et al. 2012; Schenker et al. 2012; Shibuya et al. 2012; Finkelstein et al. 2013; Tilvi et al. 2014; Huang et al. 2016; Schmidt et al. 2016; Song et al. 2016; Hoag et al. 2017, 2018; Larson et al. 2018; Pentericci et al. 2018). The likely reason for this is the increased neutral fraction at $z \gtrsim 6$ combined with the faintness of the sources (e.g., Treu et al. 2013; Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014). Interestingly, a number of recent works have reported spectroscopic confirmation for bright ($H \sim 25$ mag) LBGs at the epoch of the reionization from Ly$\alpha$ detection (Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016; Stark 2016). These observations further suggested that reionization could have happened in a patchy form, rather than homogeneously, and inspired confidence in our ability to reliably select bright sources to the highest possible redshifts.

Perhaps surprisingly, observational progress on the very bright end has been relatively slow. Covering wide areas with HST is very inefficient due to the extremely low surface densities of the brightest $z > 8$ galaxies. Some progress has come from pure parallel imaging surveys such as BORG/HIPPIES (Trenti et al. 2011; Yan et al. 2011), from targeted follow-up over the full Cosmic Assembly Near-Infrared Deep Extragalactic Survey (CANDELS) area (Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016; Stark 2016) and from the RELICS program (Salmon et al. 2017), which builds on the strong-lensing strategy of the Hubble Frontier Field (HFF) and CLASH surveys. Combined together, these wider-area, shallow surveys still only cover <1300 arcmin$^2$ and provided only $\lesssim 5$ candidates at $z \gtrsim 8$ brighter than $M_{UV} \lesssim -22.0$ (Bernard et al. 2016; Calvi et al. 2016; Livermore et al. 2018; Morishita et al. 2018).

An alternative approach consists in leveraging the ongoing wide-field ground-based surveys such as COSMOS/UltraVISTA and UKIDSS/UDS, which benefit from deep ($\sim 26$ mag) wide wavelength coverage (0.3–5 µm—e.g., Bowler et al. 2012, 2014, 2015, 2017; Stefanon et al. 2017b).

Here we report the full analysis and the results of the search for ultrabright $H \sim 24-26$ mag galaxy candidates at $z \gtrsim 8$ from the COSMOS/UltraVISTA program. This search takes advantage of the deepest-available ground-based optical+near-infrared observations, in particular the third data release (DR3) of UltraVISTA, which provides $\sim 1.4$ mag deeper data in Y, J, H, and $K_s$ compared with DR1 (McCracken et al. 2012). Our study also takes advantage of deep Spitzer/IRAC (Fazio et al. 2004) observations from the Spitzer Large Area Survey with Hyper Suprime-Cam (SPLASH, PI: Capak) and the Spitzer Matching survey of the UltraVISTA ultradeep Stripes (SMUVS, PI: Caputi–Caputi et al. 2017; Ashby et al. 2018) programs. The increased depth and the inclusion of Spitzer/IRAC data, probing the rest-frame optical, now makes it possible to access the galaxy population at $z \gtrsim 8$ through reliable sample selections.

In Stefanon et al. (2017b) we already presented five candidate bright $z \gtrsim 8$ LBGs initially identified in this search. Specifically, in that work we focused on the analysis of those sources with recent HST/WFC3 imaging from one of our programs, and showed that the new HST observations strengthened the available photometric constraints placing them at $z \sim 8$. The purpose of the present work is to present the parent sample from which those five objects were selected.

This paper is organized as follows. The observations are summarized in Section 2, while in Section 3 we describe how we performed the photometry. The source selection is detailed in Section 4. The sample is presented in Section 5 and it is characterized in Section 6. We present our conclusions in Section 7. Throughout, we adopt $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$. Magnitudes are given in the AB system from Oke & Gunn (1983), and we adopt a Chabrier (2003) initial mass function (IMF).

## 2. Observational Data

Our analysis is based on ultradeep near-infrared imaging over the COSMOS field (Scoville et al. 2007) from the DR3 of UltraVISTA. UltraVISTA provides imaging that covers 1.6 square degrees (McCracken et al. 2012) in the Y, J, H, and $K_s$ filters to $\sim 24-25$ mag (AB, $5\sigma$), with DR3 achieving fainter limits over 0.8 square degrees in 4 ultradeep stripes. The DR3 contains all data taken between 2009 December and 2014 July and reaches $Y = 25.4$, $J = 25.4$, $H = 25.1$, and $K = 24.8$ mag (AB, $5\sigma$ in $1''^2$-diameter apertures). The nominal depth we measure in the Y, J, H, and $K_s$ bands for the UltraVISTA DR3 release is $\sim 0.2$ mag, $\sim 0.6$ mag, $\sim 0.8$ mag, and $\sim 0.2$ mag, respectively, deeper than in the UltraVISTA DR2 release.

The optical data consist of CFHT/Megacam in *g*, *r*, *i*, *y*, and *z* (Erben et al. 2009; Hildebrandt et al. 2009 from the Canada-France-Hawaii Legacy Survey (CFHTLS)), Subaru/Suprime-Cam $B_j$, $V_j$, $g^+$, $r^+$, $i^+$, and $z^+$-imaging (Taniguchi et al. 2007), and Subaru Hyper Suprime-Cam (HSC) *g*, *r*, *i*, *z*, and *y* (Aihara et al. 2018a, 2018b).

For this work, we used full-depth Spitzer/IRAC 3.6 and 4.5 µm mosaics we built, combining observations from all available programs: S-COSMOS (Sanders et al. 2007), the Spitzer Extended Deep Survey (Ashby et al. 2013), the Spitzer-Cosmic Assembly Near-Infrared Deep Extragalactic Survey (S-CANDELS, Ashby et al. 2015), SPLASH (PI: Capak), and SMUVS (Caputi et al. 2017; Ashby et al. 2018). Compared with the original S-COSMOS IRAC data, SPLASH provides a large improvement in depth over nearly the whole UltraVISTA
<table><caption>**Table 1**
Photometric Depths of the Adopted Ground-based and *Spitzer*/IRAC Data Sets, and Corresponding Average Aperture Corrections</caption><thead><tr><th>Filter Name</th><th>Aperture Correction<sup>a</sup></th><th>Depth 5σ<sup>b</sup></th></tr></thead><tbody><tr><td>CFHTLS <i>u</i><sup>*</sup></td><td>2.2</td><td>26.7</td></tr><tr><td>SSC <i>B</i></td><td>1.7</td><td>27.4</td></tr><tr><td>HSC <i>g</i><sup>c</sup></td><td>2.1</td><td>26.7</td></tr><tr><td>CFHTTLS <i>g</i></td><td>2.1</td><td>26.8</td></tr><tr><td>SSC <i>V</i></td><td>2.1</td><td>26.4</td></tr><tr><td>HSC <i>r</i><sup>c</sup></td><td>1.7</td><td>26.8</td></tr><tr><td>CFHTTLS <i>r</i></td><td>2.0</td><td>26.4</td></tr><tr><td>SSC <i>r</i><sup>+</sup></td><td>2.0</td><td>26.6</td></tr><tr><td>SSC <i>i</i><sup>+</sup></td><td>1.9</td><td>26.2</td></tr><tr><td>CFHTTLS <i>y</i></td><td>1.9</td><td>26.1</td></tr><tr><td>CFHTTLS <i>i</i></td><td>1.9</td><td>26.0</td></tr><tr><td>HSC <i>i</i><sup>c</sup></td><td>1.8</td><td>26.3</td></tr><tr><td>CFHTTLS <i>z</i></td><td>2.0</td><td>25.2</td></tr><tr><td>HSC <i>z</i><sup>c</sup></td><td>1.7</td><td>25.9</td></tr><tr><td>SSC <i>z</i><sup>+</sup></td><td>2.2</td><td>25.0</td></tr><tr><td>HSC <i>y</i><sup>c</sup></td><td>2.1</td><td>24.9</td></tr><tr><td>UVISTA <i>Y</i></td><td>2.5</td><td>25.4/24.5</td></tr><tr><td>UVISTA <i>J</i></td><td>2.3</td><td>25.4/24.4</td></tr><tr><td>UVISTA <i>H</i></td><td>2.2</td><td>25.1/24.1</td></tr><tr><td>UVISTA <i>K</i><sub>s</sub></td><td>2.1</td><td>24.8/23.7</td></tr><tr><td>IRAC 3.6 µm</td><td>2.7<sup>d</sup></td><td>25.4/24.9/24.5</td></tr><tr><td>IRAC 4.5 µm</td><td>2.7<sup>d</sup></td><td>25.3/24.7/24.3</td></tr><tr><td>IRAC 5.8 µm</td><td>3.4<sup>d</sup></td><td>20.8</td></tr><tr><td>IRAC 8.0 µm</td><td>4.1<sup>d</sup></td><td>20.6</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Average multiplicative factors applied to estimate total fluxes.

<sup>b</sup> Average depth over the full field corresponding to 5σ flux dispersions in empty apertures of 1′′2 diameter corrected to total using the average aperture correction. The two depths for UltraVISTA correspond to the ultradeep and deep stripes, respectively; the three depths for the *Spitzer*/IRAC 3.6 µm and 4.5 µm bands correspond to the regions with SMUVS+SCOSMOS+SPLASH coverage (approximately overlapping with the ultradeep stripes) and SPLASH +SCOSMOS only (≈deep stripes).

<sup>c</sup> The Hyper Suprime-Cam data were not available during the initial selection of the sample; we included them in our subsequent analysis, applying the same methods adopted for the rest of the ground and *Spitzer*/IRAC mosaics.

<sup>d</sup> Aperture corrections for IRAC bands refer to the 1′′8 diameter.

area, covering the central 1.2 square degree COSMOS field to 25.5 mag (AB) at 3.6 and 4.5 µm. Spectral energy distributions (SEDs) and S-CANDELS cover smaller areas to even deeper limits, while SMUVS pushes deeper over the ultradeep UltraVISTA stripes.

Finally, we also included measurements in the IRAC 5.8 µm and 8.0 µm bands from the S-COSMOS program. Even though the coverage in these bands is rather shallow (~20.7 mag, 5σ in 1′′8-diameter aperture), detections in these two bands can be useful to discriminate high-redshift sources from lower-redshift interlopers. We discuss this for our sample at the end of Section 5.2.

A summary of all the deep, wide-area data sets along with 5σ depths is provided in Table 1, while in Figure 1 we present the coverage of the different data sets.

### 3. Photometry

Source catalogs were constructed using SExTRACTOR v2.19.5 (Bertin & Arnouts 1996), run in dual image mode, with source detection performed on the square root of a χ²

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Depth and layout of observations relevant to our current search for $z \sim 8$ and $z \sim 9$ galaxies over the UltraVISTA field. The gray shaded image represents the UltraVISTA DR3 exposure time map (deeper exposure for darker regions). The colored curves mark the coverage from CFHT Legacy Survey (magenta), ultra-deep HSC (green), and the deep *Spitzer*/IRAC observations from the SPLASH program (red). Even deeper *Spitzer*/IRAC observations are available over the deep stripes from the SMUVS program. The yellow rectangle in the center demarcates the region with observations from the CANDELS program. The blue-shaded image corresponds to the COSMOS/DASH coverage map (darker regions indicate deeper coverage). The orange stars mark the position of bright candidate $z \sim 8$ galaxies we have discovered in our search.</figcaption></figure>

image (Szalay et al. 1999) built from the combination of the UltraVISTA *J*, *H*, and $K_s$ images.

The first selection was performed adopting ground-based observations only. Images were first convolved to the *J*-band point-spread function (PSF) and carefully registered against the detection image (mean rms ~0.05). Initial color measurements were made in small Kron (1980)-like apertures (SExtractor AUTO and Kron factor 1.2) with typical radius $r_{color} \sim 0.35-0.50$.

Successively, we refined our selection of $z \sim 8$ and $z \sim 9$ candidate galaxies using color measurements made in fixed 1′′2-diameter apertures. For this step, fluxes from sources and their nearby neighbors (12′′0 × 12′′0 region) are carefully modeled; aperture photometry is then performed after subtracting the neighbors using MOPHONGO (Labbé et al. 2006, 2010a, 2010b, 2013, 2015). Our careful modeling of the light from neighboring sources improves the overall robustness of our final candidate list to source confusion. Total magnitudes are derived by correcting the fluxes measured in 1′′2-diameter apertures for the light lying outside this aperture. The relevant correction factor is estimated on a source-by-source basis based on the spatial profile of each source and the relevant PSF-correction kernel. Average PSF corrections for each band are listed in Table 1.

Photometry on the *Spitzer*/IRAC observations is more involved due to the much lower-resolution FWHM = 1′′7 compared to the ground-based data (FWHM = 0.7). The lower resolution results in source blending, where light from foreground sources contaminates measurements of the sources of interest. Photometry of the IRAC bands was therefore performed with MOPHONGO, adopting 1′′8 apertures. Similarly to the optical bands, IRAC fluxes were corrected to total for missing light outside the aperture using the model profile for
the individual sources. The procedure for IRAC photometry employed here is very similar to those of other studies (e.g., Galametz et al. 2013; Guo et al. 2013; Skelton et al. 2014; Nayyeri et al. 2017; Stefanon et al. 2017a).

Following Stefanon et al. (2017b), the uncertainties associated with the flux densities were estimated from the standard deviation of the flux density measurements in 1"2-diameter empty apertures, multiplied by the corresponding aperture correction.

## 4. Sample Selection

We require sources to be detected at >5σ significance in the J, H, Ks, [3.6], and [4.5] images after coadding their signal-to-noise ratios (S/Ns) in quadrature and in those bands with a positive flux density estimate, and we limit our selection to sources brighter than H ~ 25.8 mag. The combined Ultra-VISTA and IRAC detection and S/N requirements exclude spurious sources due to noise, detector artifacts, and diffraction features.

We identified candidate z ~ 8 and z ~ 9 LBGs using a combination of Lyman-break criteria and photometric redshift selections. While photometric redshifts are a great tool in a number of cases, their quality is a direct consequence of the adopted set of template models. It is not uncommon, for instance, when running photometric redshift codes to obtain solutions at z ≳ 6 represented by red, dusty SEDs. Given our current limited knowledge on the physical properties of high-redshift galaxies, the existence of such objects, although unlikely, is still possible. However, their red colors would make the assessment of their nature very difficult with the available data, being unable to effectively exclude (more likely) low-redshift solutions. The LBG cuts we applied are strict enough to exclude sources with red, power-law-like SEDs, therefore aiming at selecting the most robust sample of star-forming galaxies consistent with at most a small amount of dust attenuation. Furthermore, because the process we applied to measure flux densities heavily relies on MOPHONGO, it would have required an unfeasible amount of time, running it on 24 bands for the full set of sources detected on the χ² image (~1 million sources). For these reasons, we started from a subsample selected with Lyman-break cuts, and consolidated the selection, applying a photometric redshift analysis. The full procedure is detailed below.

We construct a preliminary catalog of candidate z ~ 8 and z ~ 9 galaxies using those sources that show an apparent Lyman-break due to absorption of UV photons by neutral hydrogen in the intergalactic medium (IGM) blueward of the redshifted Lyα line. At z > 7.1, the break results in a significantly lower Y-band flux density for candidates, while at z > 8.7 it reduces the J-band flux densities. Because of this we applied two distinct criteria to select either z ~ 8 or z ~ 9 candidate LBGs. Specifically, for the z ~ 8 sample we applied the following criterion:

$$Y - (J + H)/2 > 0.75 \qquad (1)$$

while for the z ~ 9 sample we required that:

$$J - H > 0.8. \qquad (2)$$

In case of a non-detection, the Y- or J-band flux in these relations was replaced by the equivalent 1σ upper limit.

These cuts do not exclusively select z > 7 galaxies, but also accept some dust-reddened low-redshift galaxies. However,

such sources would show a very red continuum and red colors redward of the J-band or H-bands. Therefore, to reject this class of galaxies we also imposed to each one of the sample selected with Equations (1) and (2) the requirement of a blue continuum redward of the break:

$$(H - K < 0.7) \wedge ((K - [3.6] < 1.75) \lor (H - [3.6] < 1.75)) \qquad (3)$$

where ∧ denotes the logical AND operator, and ∨ denotes the logical OR operator. These limits are valuable for excluding a small number of very red sources from our selection. Nevertheless, it is worth emphasizing that our final sample of z > 7 bright galaxies shows little dependence on the specific limits chosen here. This initial selection resulted in 2234 candidates (out of ~8 × 10⁵ detected sources), 2015 Y-dropouts, and 183 J-dropouts.

We further cleaned our sample from low-redshift sources and Galactic stars by imposing χ²opt < 4. The χ²opt is defined as χ²opt = ΣᵢSGN(fᵢ)(fᵢ/σᵢ)² (Bouwens et al. 2011), where fᵢ is the flux in any optical band i with uncertainty σᵢ, and SGN(fᵢ) is +1 if fᵢ > 0 and -1 if fᵢ < 0. The χ²opt is calculated in both 1"2-diameter apertures and in the scaled elliptical apertures. χ²opt is effective in excluding z = 1–3 low-redshift star-forming galaxies where the Lyman-break color selection is satisfied by strong line emission contributing to one of the broad bands (e.g., Atek et al. 2011; van der Wel et al. 2011). We also constructed full-depth pseudo r-, i-, and z-band mosaics, combining the relevant observations from the CFHTLS, HSC, and Subaru SuprimeCam (SSC) data sets, and excluded sources with a 2σ detection in either individual ground-based imaging bands or in one of the three full-depth optical mosaics, as potentially corresponding to lower-redshift contaminants. This step left 901 candidate LBGs in our sample (791 Y-dropouts and 110 J-dropouts).

Subsequently, we determined the redshift probability distribution P(z). For this we used the EAZY program (Brammer et al. 2008), which fits non-negative linear combination of galaxy spectral templates to the observed SED, assuming a flat prior on redshifts. We complemented the standard EAZY_v1.0 template set with templates extracted from the Binary Population and Spectral Synthesis code (BPASS—Eldridge et al. 2017) v1.1 for subsolar metallicity (Z = 0.2Z⊙), which include nebular emission from CLOUDY. Specifically, we adopted templates with equivalent widths EW(Hα) ~ 1000–5000 Å, as these extreme EW reproduce the observed [3.6]–[4.5] colors for many spectroscopically confirmed z ~ 7–9 galaxies (Ono et al. 2012; Finkelstein et al. 2013; Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016; Stark 2016). Driven by current observational results (e.g., Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016), we blanketed the Lyα line from those templates with EW(Lyα) ≳ 40 Å. Finally, we added templates of 2 Gyr old, passively evolving systems from Bruzual & Charlot (2003), with Calzetti et al. (2000) extinction in the range A$_{V}$ = 0–8 mag to test the robustness of our selected candidates against being lower-redshift interlopers highly attenuated by dust. We imposed an additional constraint: that the integrated probability beyond z = 6 be >50%. The use of a redshift likelihood distribution P(z) is very effective in rejecting faint low-redshift galaxies with a strong Balmer/4000 Å break and fairly blue colors redward of the break. After this step,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Expected redshift distributions (normalized to unit area) for Y- and J-dropouts recovered from our Monte Carlo simulation (Section 6.5) for the UltraVISTA sample (blue and magenta solid curves, respectively). For comparison, we also show the expected redshift distributions for the Y-dropout LBGs samples from Bouwens et al. (2015). The two $z \sim 8$ distributions largely overlap, supporting the combination of the UltraVISTA with Bouwens et al. (2015) LFs.</figcaption></figure>

the sample resulted in 49 candidates (44 Y-dropouts and 5 J-dropouts).

In Figure 2 we present the expected redshift distribution of the Y-(i.e., $z \sim 8$) and J- ($z \sim 9$) dropout selections obtained from our Monte Carlo simulations described in Section 6.5. The Y-dropout selection over UltraVISTA peaks at $z \sim 8.2$, but the wings of the $z_{phot}$ extend into $z \sim 9$. On the other side, the distribution of $z_{phot}$ from the J-dropout selection presents a wing at lower redshifts, reaching $z \sim 7.5$ introduced by the lack of continuity in the coverage of wavelengths between the Y- and J-bands from the atmospheric absorption.

All the 49 candidates showed compact morphologies. However, the relatively low S/N and coarser spatial resolution of the ground-based data make the distinction between a point source (indicative of a low-mass star nature) and an extended object challenging. Therefore, to further exclude contamination by the coolest low-mass stars we used EAZY to fit all candidates with stellar templates from the SpecX prism library (Burgasser 2014) and exclude any that are significantly better fit ($\Delta\chi^2 > 1$) by stellar SED models. The approach we utilized is identical to the SED-fitting approach recently employed by Bouwens et al. (2015) for excluding low-mass stars from the CANDELS fields. Through this step we excluded 30 sources as likely brown-dwarf candidates. The sample surviving this selection included 17 candidate Y-dropout LBGs and 2 candidates J-dropouts.

The IRAC flux densities are particularly crucial for our work, because of the dependence of the [3.6]–[4.5] color on redshift, and because for $z \gtrsim 8$ the 3.6 µm and 4.5 µm bands probe the rest-frame optical redward of the Balmer break, thus providing information of the age and stellar mass of the sources. For these reasons, we visually inspected the image stamps containing the original IRAC science frame subtracted of the model sources (hereafter *residual images*). Residual images showed generally clean subtractions, with the exception of two sources (UVISTA-Y7 and UVISTA-Y9). Because the photometric redshifts for these two sources obtained after excluding the IRAC bands still indicated a $z \sim 8$ solution, we opted for including the two sources when estimating the LF (see Section 6.5), but we excluded them from physical

parameter considerations due to likely suffering from systematics (Sections 6.2–6.4).

Finally, we excluded one Y-dropout source that, even though it satisfied all the previous criteria, showed a $2.2\sigma$ detection on the image built stacking all the optical data.

When considered together, our selection criteria resulted in very low expected contamination rates. The nominal contamination rate just summing over the redshift likelihood distribution for the $z \sim 8$ sample is $\sim 5\%$, based on the assumption our SED templates span the range of colors for the low-$z$ interlopers. This percentage should just be considered indicative; it does not account for $z < 6$ sources scattering into our selection due to the impact of noise. We will conduct such a quantification in Section 5.4.

In addition to minimizing the impact of contamination in our $z \sim 8$ selection, the present selection criteria also likely exclude some bona fide $z \sim 8$ galaxies and thus introduce some incompleteness into our $z \gtrsim 8$ samples. We cope with this incompleteness using selection volume simulations in Section 6.5.

## 5. Results

The above selection criteria resulted in a total of 18 LBG candidates over the UltraVISTA field. Specifically, we identified 16 Y-band dropouts (likely $z \sim 8$ candidate LBGs) and 2 J-band dropouts (likely $z \sim 9$ candidate LBGs). These candidates span a range of $H \sim 24.0–26.0$ mag and constitute the most luminous $z \sim 8$ galaxy candidates known to date, 0.5–1.0 mag brighter than the galaxies recently confirmed through spectroscopy (Oesch et al. 2015b; Zitrin et al. 2015; Roberts-Borsani et al. 2016).

Stefanon et al. (2017b) already presented five of them: three Y-band dropouts (namely UVISTA-Y1, UVISTA-Y5, and UVISTA-Y6) and the two J-band dropouts (UVISTA-J1 and UVISTA-J2) that we had followed up with HST/WFC3 imaging in the F098W, F125W, and F160W bands. That analysis further supported the conclusion that the three Y-band dropouts are $z \gtrsim 8$ LBGs, and showed that the two J-band dropout candidates were low-redshift interlopers. In the next sections we present the full sample from which those five sources were extracted. For completeness, we also reexamined the three sources analyzed in Stefanon et al. (2017b) (UVISTA-Y1, UVISTA-Y5, and UVISTA-Y6), excluding the flux density measurements in the HST/WFC3 bands, and conclude that they are probable $z \gtrsim 8$ candidates. We refer the reader to Stefanon et al. (2017b) for full details on their analysis, including the HST flux densities. Nonetheless, high-resolution imaging from HST is key in ascertaining the nature of these sources, as we discuss in the next section.

### 5.1. *High-resolution Imaging from HST*

In an effort to further ascertain the nature of the $z \sim 8$ LBG sample considered in this work, we also inspected the recent Drift And SShift mosaic (DASH—Momcheva et al. 2016; Mowla et al. 2019) at the nominal locations of the selected candidate bright LBGs. This mosaic covers $\sim 0.7$ sq. deg of sky in the WFC3/F160W band to a depth of $\sim 25.1$ mag (0.′′3 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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Image stamps (5′′0 side) of those sources with coverage in the WFC3/F160W DASH mosaic (Momcheva et al. 2016; Mowla et al. 2018), centered at the nominal location of each object. To improve contrast, each cutout has been smoothed with a 0.′′2 Gaussian filter.</figcaption></figure>

candidate LBGs was performed on ground-based data (seeing FWHM ~ 0.′′7), the finer spacial resolution of HST/WFC3 (PSF FWHM ~ 0.′′2) 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.

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 $\gtrsim 4\sigma$), 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*<sub>s</sub> 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 ~0.′′5. 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 \gtrsim 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 \sim 8$ source initially identified on the UltraVISTA images resulted in the three objects being at $z \gtrsim 8$.

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

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Image stamp (5′′0 × 5′′0, smoothed with a Gaussian of 0.′′1 FWHM) in the WFC3/F160W band extracted from the DASH mosaic (Momcheva et al. 2016; Mowla et al. 2018) centered at the position of UVISTA-Y3. Individual components are indicated by the blue labels. The red curve corresponds to the contour of the stacked *J*, *H*, and *K*<sub>s</sub> data.</figcaption></figure>

<table><caption>**Table 2**
Candidate $z \sim 8$ LBGs with HST/WFC3 F160W Coverage</caption><thead><tr><th>ID</th><th>PID</th><th>PI</th><th>Depth (mag)</th></tr></thead><tbody><tr><td>UVISTA-Y1</td><td>14895</td><td>R. Bouwens</td><td>24.7</td></tr><tr><td>UVISTA-Y2</td><td>14114</td><td>P. van Dokkum</td><td>24.9</td></tr><tr><td>UVISTA-Y3<sup>a</sup></td><td>13868</td><td>D. Kocevski</td><td>26.5</td></tr><tr><td>UVISTA-Y4</td><td>14114</td><td>P. van Dokkum</td><td>24.9</td></tr><tr><td>UVISTA-Y5</td><td>14895</td><td>R. Bouwens</td><td>24.9</td></tr><tr><td>UVISTA-Y6</td><td>14895</td><td>R. Bouwens</td><td>25.0</td></tr><tr><td>UVISTA-Y7</td><td>14114</td><td>P. van Dokkum</td><td>24.9</td></tr><tr><td>UVISTA-Y8</td><td>13641</td><td>P. Capak</td><td>25.7</td></tr><tr><td>UVISTA-Y9</td><td>14114</td><td>P. van Dokkum</td><td>24.8</td></tr><tr><td>UVISTA-Y10</td><td>14114</td><td>P. van Dokkum</td><td>24.7</td></tr><tr><td>UVISTA-Y11</td><td>12440</td><td>S. Faber</td><td>26.6</td></tr><tr><td>UVISTA-Y13</td><td>14114</td><td>P. van Dokkum</td><td>24.9</td></tr><tr><td>UVISTA-Y14</td><td>14114</td><td>P. van Dokkum</td><td>24.8</td></tr><tr><td>UVISTA-Y16</td><td>14114</td><td>P. van Dokkum</td><td>24.7</td></tr></tbody></table>

**Note.** The limiting magnitudes refer to $5\sigma$ fluxes in apertures of 0.′′6 diameter corrected to total using the growth curve of point sources.

<sup>a</sup> HST/WFC3 imaging suggests this source is potentially multiple. See Section 5.1 for details.

are $r_e \lesssim 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 \sim 8$, a separation of 0.′′5 corresponds to ~2.5 kpc, i.e., $\gtrsim 2.5 \times$ 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 \sim 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.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Image stamps ($5''0 	imes 5''0$) for UVISTA-Y11 in HST bands from the CANDELS program, as labeled at the top-left corner of each panel. No evidence for flux at the nominal location of the source is seen blueward of the 1.2 µm band, consistent with what is seen in the ground-based observations.</figcaption></figure>

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 \sim 8$ candidates over the $\sim 0.8$ deg$^2$ of the UltraVISTA ultradeep stripes, we would expect to find only $\sim 1$ candidate over the $\sim 190$ arcmin$^2$ CANDELS COSMOS field. Indeed, only one $z \sim 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_{606}$ $I_{814}$, $J_{125}$, $JH_{140}$, and $H_{160}$. The $V_{606}$ mosaic shows a close low-$z$ neighbor just $\sim 0''7$ 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.$^{14}$ We do not detect flux at $>1\sigma$ in the $V_{606}$ and $I_{814}$ bands, increasing our confidence in its high-$z$ nature.

Finally, we inspected the ACS $I_{814}$-band mosaic of the COSMOS program (Scoville et al. 2007, $\sim 26.5$ mag in $0''6$ aperture diameter, $5\sigma$). 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 $\sim 1''0$ 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 \sim 8$ LBGs stresses the importance of deep ($\gtrsim 1$ orbit) high-resolution multi-band follow-up to further assess the nature of the remarkable LBG candidates identified in the present work.

## 5.2. Sample of $z \sim 8$ Candidates

Figure 6 presents the image stamps of all the candidate $z \sim 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\sigma$, adding in quadrature the detection significance in the J, H, and Ks bands.

The observed SEDs of the galaxy candidates are presented in Figure 7, along with the EAZY best-fit templates at $z \sim 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 \sim 8$ solution, consistent with our initial selection.

Four of our 16 $z \sim 8$ candidates (or $\sim 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 $\sim 37\%$ of the area, we find no evidence for a higher surface density of $z \sim 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 \sim 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 \gtrsim 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 \lesssim p(z > 7) \lesssim 0.95$.$^{15}$ 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 \sim 8$ selection (UVISTA-Y3 and UVISTA-Y13) are formally detected at $>1\sigma$ in the combined 5.8 µm and 8.0 µm observations, with nominal brightnesses of $\sim 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 ($\sim 25$ mag, i.e., an approximately flat $f_\nu$ SED) at $\sim 7$µm, simple noise statistics predict $4 \pm 2$ sources to be detected at $>1\sigma$. We therefore conclude that the

<footnote>$^{14}$ Omitting the neighboring source leads to flux densities systematically overestimated by $\sim 30\%$.</footnote>

<footnote>$^{15}$ When considered as a single object, UVISTA-Y3 has a $p(z > 7) = 0.99$, suggesting a fairly robust redshift for this source as well.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Stacked ground-based optical, near-infrared, and *Spitzer*/IRAC image stamps for our bright candidate *z* ~ 8 galaxies selected over COSMOS/UltraVISTA. Each image stamp is 10''0 × 10''0 in size and it is shown in inverted grayscale. Neighbor-subtraction was applied to the IRAC data.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. (Continued.)</figcaption></figure>

>1σ formal detection of two $z \sim 8$ candidates in our selection is not a concern.

### 5.3. Sample of $z \sim 9$ Candidates

The selection criteria expressed by Equations (2) and (3) are designed to select $z \gtrsim 9$ LBG candidates. Indeed our initial analysis identified two exceptionally bright ($m_H \sim 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 \sim 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 \sim 8$, in the rest of this work we consider the 16 Y-band dropouts as our fiducial sample of galaxies at $z \sim 8$; specifically, we include in the $z \sim 8$ sample those Y-dropouts with nominal $z_{phot} \sim 9$ (see Figure 2). However, in Section 6.5 we also consider the contribution of those sources with $z_{phot} \sim 9$ to the $z \sim 9$ LF. We refer the reader to our discussion in Section 6.5 for full details.

### 5.4. Expected Contamination in Our Bright $z \sim 8$ Samples

One potentially important source of contamination for our current $z \sim 8$ and $z \sim 9$ samples occurs through the impact of
<table><caption>**Table 3**
Sample of Candidate *z* ~ 8 LBGs</caption><thead><tr><th>ID</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>m<sub>H</sub><sup>a</sup><br>(mag)</th><th>Y - J<sup>b</sup><br>(mag)</th><th>[3.6]-[4.5]<sup>b</sup><br>(mag)</th><th>z<sub>phot</sub><sup>c</sup></th></tr></thead><tbody><tr><td>UVISTA-Y1<sup>d,i</sup></td><td>09: 57: 47.900</td><td>+02: 20: 43.66</td><td>24.8 ± 0.1</td><td>>2.1</td><td>0.4 ± 0.2</td><td>8.53<sup>+0.53</sup><sub>-0.62</sub></td></tr><tr><td>UVISTA-Y2<sup>i</sup></td><td>10: 02: 12.558</td><td>+02: 30: 45.71</td><td>24.8 ± 0.2</td><td>>2.2</td><td>0.5 ± 0.1</td><td>8.21<sup>+0.50</sup><sub>-0.49</sub></td></tr><tr><td>UVISTA-Y3a<sup>c</sup></td><td>10: 00: 32.324</td><td>+01: 44: 30.86</td><td>25.5 ± 0.3</td><td>>0.9</td><td>0.6 ± 0.5</td><td>8.68<sup>+0.93</sup><sub>-1.21</sub></td></tr><tr><td>UVISTA-Y3b<sup>c</sup></td><td>10: 00: 32.317</td><td>+01: 44: 31.48</td><td>26.1 ± 0.5</td><td>>0.9</td><td>&lt;0.8<sup>f</sup></td><td>8.90<sup>+1.24</sup><sub>-1.18</sub></td></tr><tr><td>UVISTA-Y3c<sup>e</sup></td><td>10: 00: 32.350</td><td>+01: 44: 31.73</td><td>26.0 ± 0.5</td><td>>-0.5</td><td>0.7 ± 0.5</td><td>9.29<sup>+1.58</sup><sub>-2.10</sub></td></tr><tr><td>UVISTA-Y4<sup>i</sup></td><td>10: 00: 58.485</td><td>+01: 49: 55.96</td><td>24.9 ± 0.2</td><td>1.0 ± 0.4</td><td>0.1 ± 0.2</td><td>7.42<sup>+0.19</sup><sub>-0.20</sub></td></tr><tr><td>UVISTA-Y5<sup>d,i</sup></td><td>10: 00: 31.886</td><td>+01: 57: 50.23</td><td>24.9 ± 0.2</td><td>>1.3</td><td>0.8 ± 0.3</td><td>8.60<sup>+0.58</sup><sub>-0.65</sub></td></tr><tr><td>UVISTA-Y6<sup>d</sup></td><td>10: 00: 12.506</td><td>+02: 03: 00.50</td><td>25.3 ± 0.3</td><td>>1.5</td><td>0.3 ± 0.4</td><td>8.32<sup>+0.66</sup><sub>-0.92</sub></td></tr><tr><td>UVISTA-Y7<sup>i</sup></td><td>09: 59: 02.566</td><td>+02: 38: 06.05</td><td>25.5 ± 0.4</td><td>>1.3</td><td>...</td><td>8.47<sup>+0.72</sup><sub>-0.73</sub></td></tr><tr><td>UVISTA-Y8<sup>i</sup></td><td>10: 00: 47.544</td><td>+02: 34: 04.84</td><td>25.4 ± 0.3</td><td>>1.4</td><td>1.0 ± 0.8</td><td>8.34<sup>+0.60</sup><sub>-0.58</sub></td></tr><tr><td>UVISTA-Y9</td><td>09: 59: 09.621</td><td>+02: 45: 09.68</td><td>25.4 ± 0.3</td><td>0.8 ± 0.7</td><td>...</td><td>7.69<sup>+0.99</sup><sub>-0.71</sub></td></tr><tr><td>UVISTA-Y10<sup>i</sup></td><td>10: 01: 47.495</td><td>+02: 10: 15.37</td><td>25.3 ± 0.3</td><td>>1.6</td><td>0.9 ± 0.7</td><td>8.25<sup>+0.61</sup><sub>-0.60</sub></td></tr><tr><td>UVISTA-Y11<sup>i</sup></td><td>10: 00: 19.607</td><td>+02: 14: 13.15</td><td>25.2 ± 0.3</td><td>>1.4</td><td>0.8 ± 0.4</td><td>8.64<sup>+0.66</sup><sub>-0.72</sub></td></tr><tr><td>UVISTA-Y12<sup>i</sup></td><td>10: 00: 15.975</td><td>+02: 43: 32.96</td><td>25.6 ± 0.4</td><td>>1.2</td><td>0.2 ± 0.8<sup>f</sup></td><td>8.70<sup>+0.61</sup><sub>-0.74</sub></td></tr><tr><td>UVISTA-Y13</td><td>09: 58: 45.561</td><td>+01: 53: 41.79</td><td>25.8 ± 0.4</td><td>>1.1</td><td>0.8 ± 0.7</td><td>8.54<sup>+0.79</sup><sub>-1.18</sub></td></tr><tr><td>UVISTA-Y14</td><td>10: 00: 12.568</td><td>+01: 54: 28.50</td><td>25.6 ± 0.4</td><td>>1.1</td><td>0.1 ± 0.6</td><td>7.55<sup>+1.71</sup><sub>-2.68</sub></td></tr><tr><td>UVISTA-Y15</td><td>09: 57: 35.795</td><td>+02: 11: 57.81</td><td>25.6 ± 0.4</td><td>1.1 ± 0.9</td><td>&lt;-0.5<sup>f,g</sup></td><td>7.64<sup>+1.13</sup><sub>-1.13</sub></td></tr><tr><td>UVISTA-Y16<sup>i</sup></td><td>10: 01: 56.333</td><td>+02: 34: 16.25</td><td>25.3 ± 0.3</td><td>1.2 ± 0.7</td><td>0.6 ± 0.4</td><td>7.90<sup>+0.74</sup><sub>-0.57</sub></td></tr></tbody></table>

**Notes.** Measurements for the ground-based bands are 1''2 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 1''8 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.

<sup>a</sup> H-band magnitude and associated $1\sigma$ uncertainty estimated from the UltraVISTA DR3 mosaic.

<sup>b</sup> Upper/lower limits to be intended as $1\sigma$. In computing these colors, we replaced negative fluxes with their corresponding $1\sigma$ uncertainty. See Tables 7–9 for a complete listing of flux densities in all bands.

<sup>c</sup> Photometric redshift and 68% confidence interval of the best-fitting template from EAZY.

<sup>d</sup> These 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.

<sup>e</sup> These 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 \pm 0.1$ mag; $Y - J = 1.1 \pm 0.4$ mag; [3.6]–[4.5] = 0.3 ± 0.1 mag and $z_{phot} = 7.62^{+0.14}_{-0.28}$.

<sup>f</sup> This IRAC color is based on <2$\sigma$ flux density estimate in both bands.

<sup>g</sup> A blue [3.6]–[4.5] < 0 mag color might be indicative of a redshift $z \lesssim 7$.

<sup>h</sup> After 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.

<sup>i</sup> These sources have a probability $p(z > 7) \ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.

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 \sim 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 \le 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_{160} = 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 \sim 8$ samples. Over the ultradeep stripes where 95% of the sources in our $z \sim 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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spectral energy distributions from the observed ground-based optical, infrared, and *Spitzer*/IRAC photometry (filled red squares with error bars and black $2\sigma$ upper limits). The red arrows mark $2\sigma$ upper limits in the combined HSC, CFHTLS, and SSC *g*, *r*, and *i* bands. The solid blue curve corresponds to the best-fit SED provided by EAZY, while the gray line shows the best-fit SED when the fit is forced to a $z < 6$ solution. The corresponding redshifts are labeled in matching color, together with the total $\chi^2$. The inset plot on the upper-left corner of each panel presents the redshift probability distributions $P(z)$ for each candidate $z \sim 8$ galaxy.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7. (Continued.)</figcaption></figure>

$z \sim 8$ samples. The typical H-band magnitude of the expected contaminants ranged from $H \sim 25$ to 25.5 mag.

5.5. Possible Lensing Magnification

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 \ge 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 \sim 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 \sim 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
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Spectral energy distribution of UVISTA-Y3 when we do not deblend its photometry using the higher spatial resolution provided by COSMOS/DASH, but instead consider it as a single source. Same plotting conventions as in Figure 7. The solution is still a $z \sim 8$ LBG, consistent with our initial selection.</figcaption></figure>

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 ($\lesssim 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 \sim 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_{814}$-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 \sim 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 \gtrsim 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 \sim 0.5$), with a position angle of the LBG relative to the main axes of the two ellipses of $\sim 29^\circ.6$ and $\sim 4^\circ.5$, 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 $\sim 0.08$ and $\sim 0.04$ mag difference. Given the small contribution to the magnification estimates,

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Inferred UV luminosities and redshifts for the present sample of bright $z \sim 8$ LBGs (filled red circles). The red cross at the bottom-right corner corresponds to the median uncertainties for our sample. For comparison, we also present galaxies from the high-redshift samples compiled by Bouwens et al. (2015, 2016, CANDELS/*DF), Calvi et al. (2016), Bowler et al. (2017), Ono et al. (2018), Livermore et al. (2018), and Morishita et al. (2018). The blue line marks the evolution of the characteristic magnitude of the UV LF of Bouwens et al. (2015) up to $z = 8$ (solid blue line) and its extrapolation to $z \sim 11$ (dashed blue line). Our candidate LBGs lie at the high-luminosity end of all candidate $z \sim 8$ star-forming galaxies discovered to date, being at least $\sim 0.5-1.0$ mag brighter than the typical bright star-forming galaxy identified over the CANDELS fields.</figcaption></figure>

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 ($\sim 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 $\sim 1.4\times$, $\sim 1.16\times$, and $\sim 1.14\times$ from a $10^{10.7} M_\odot$, $z = 1.76$ galaxy (10:00:12.51, 02:02:57.3), $10^{10.6} M_\odot$, $z = 1.6$ galaxy (10:00:12.15, 02:02:59.6), and a $10^{10.3} M_\odot$, $z = 1.65$ galaxy (10:00:12.18, 02:03:00.7), respectively, that lie within $4''.9$, $3''.2$, and $5''.4$ of this source. Their velocity dispersions are estimated to be $259 \text{ km s}^{-1}$, $225 \text{ km s}^{-1}$, and $206 \text{ km s}^{-1}$, respectively.

**UVISTA-Y8:** This source is estimated to be amplified by $1.39\times$ from a $10^{10.8} M_\odot$ ($264 \text{ km s}^{-1}$), $z = 1.33$ galaxy (10:00:47.68, 02:34:08.4) that lies within $4''1$ of this source.

**UVISTA-Y9:** This source is estimated to be amplified by $1.37\times$ and $1.43\times$ by a $10^{11.0} M_\odot$ ($265 \text{ km s}^{-1}$), $z = 0.91$ galaxy (09:59:09.35, 02:45:11.8) and $10^{11.0} M_\odot$ ($268 \text{ km s}^{-1}$), $z = 0.93$ galaxy, respectively, that lie within $5''0$ and $4''6$ of the source.

**UVISTA-Y13:** This source is estimated to be amplified by $1.6\times$ by a $10^{11.15} M_\odot$ ($330 \text{ km s}^{-1}$), $z = 1.63$ galaxy (09:58:45.83, 01:53:40.6) that lies within $4''2$ 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.

## 6. Discussion

### 6.1. Bright Candidate LBGs at $z \sim 8$

In Figure 9 we present our sample of candidate $z \sim 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.

## 6.2. Rest-frame Colors of Bright $z \sim 8$ Galaxies

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 $\beta$ (where $\beta$ 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 $\beta$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$, $\beta$s derived in this way could still be affected by the Ly$\alpha$ emission line shifting into the J-band. We note, however, that observed Ly$\alpha$ equivalent widths of bright $z \sim 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 ($\beta_{phot}$) resulted in values essentially equal to those from the method we initially applied (median $\beta_{phot}-\beta_{BC03} \sim 0.1$), although with large scatter for ~30% of the sources ($\Delta\beta \gtrsim 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 $\beta$ and rest-frame $u-g$ colors for the bright $z \sim 8$ sample. The $z \sim 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, $\beta = -2.2 \pm 0.6$, is bluer but still consistent with the UV-continuum slopes found for bright $-22 < M_{UV} < -21$ galaxies at $z = 6$ ($\beta = -1.55 \pm 0.17$) and $z \sim 7$ ($\beta = -1.75 \pm 0.18$) by Bouwens et al. (2014) and suggests a continuing trend toward bluer $\beta$s at higher redshifts.

Recently, Oesch et al. (2013) analyzed the rest-frame UV and optical properties of a sample of $z \sim 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} \sim -18$ to

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Distribution of UV-continuum slopes and rest-frame $u-g$ colors for the bright $z \sim 8$ sample. The vector at the bottom-right corner shows the impact of adding a Calzetti et al. (2000) extinction of $A_V = 0.5$ mag. The scatter of points likely reflects a mixture of intrinsic variation and measurement uncertainties. There is no apparent correlation between $\beta$ and rest-frame $u-g$ as might be expected if dust were primarily responsible for the variation in both colors.</figcaption></figure>

$\sim -22$ AB. Their $J_{125} - [4.5]$ color (corresponding to approximately rest-frame $u-z$ at $z \sim 4$) shows a correlation with the UV slope $\beta$ (see e.g., their Figure 4), likely driven by dust extinction. The uniform scatter observed at $z \sim 8$, then, may suggest rapidly evolving physical mechanisms responsible for the production of dust during the ~800 Myr between the two epochs.

## 6.3. Constraints on the EWs of the [O III]+H$\beta$ Lines

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$\beta$ 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 \sim 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_{\nu}$ (i.e., a continuum [3.6]–[4.5] = 0 mag or $f_{\lambda} \propto \lambda^{-2}$), with the strongest emission lines ([O II]$_{3727}$, H$\beta$, [O III]$_{4959,5007}$, H$\alpha$, [N II]$_{6548,6583}$, [S II]$_{6716,6730}$), empirical emission lines ratios from Anders & Alvensleben (2003) for 0.2 Z$_{\odot}$ metallicity, they inferred a median [O III]+H$\beta$ 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$\alpha$
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (Left) Observed [3.6]–[4.5] colors vs. photometric redshift for our $z \sim 8$ sample (blue circles) and those from Roberts-Borsani et al. (2016, yellow squares). The predicted dependence of the [3.6]–[4.5] color on redshift is also shown for $H\alpha$ EWs of 200 Å (red), 1000 Å (purple), and 2000 Å (blue). (Right) Number of sources in our $z \sim 8$ sample (blue histogram) and that of Roberts-Borsani et al. (2016, yellow histogram) with a given [3.6]–[4.5] color. The median [3.6]–[4.5] color is 0.62 mag. On the upper horizontal axis, we present the EW([O III+$H\beta$]) corresponding to a given [3.6]–[4.5] color, assuming an intrinsic stellar continuum color of 0 mag.</figcaption></figure>

emission. So it is unclear if those results were representative of the general bright $z \sim 8$ population.

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 \sim 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\beta$ 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 \sim 8 \ L < L^* \ \text{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]$_{5007}$, depend strongly on metallicity (e.g., Inoue 2011). Considering this, we estimate that simple approximations are probably uncertain by factors of 2–3.

## 6.4. Stellar Populations of Bright $z \sim 8$ Galaxies

In this section we present our estimates of stellar population parameters for the bright $z \sim 8$ galaxies. Measurements were performed with the FAST code (Kriek et al. 2009), adopting Bruzual & Charlot (2003) models for subsolar $0.2Z_{\odot}$ 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\beta$ from the models calculated by Inoue (2011). The luminosity in $H\beta$ 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]$_{5007}$. 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$_{1600}$ absolute magnitude, the UV-continuum slope $\beta$, and the rest-frame $u-g$ color for each individual candidate bright $z \sim 8$ LBG. A summary of the physical properties is presented in Table 5.

As we already introduced in Section 4, the neighbor-cleaned IRAC 3.6 and 4.5 µm band image sections for two sources
<table><caption>**Table 4**
Main Physical Parameters for the Sample of Candidate *z* ~ 8 LBGs</caption><thead><tr><th>ID</th><th><i>M</i><sub>UV</sub><br>(mag)</th><th>UV Slope β</th><th><i>u</i> − <i>g</i><br>(mag)</th><th>log(<i>M</i><sub>*</sub>)<br>(<i>M</i><sub>☉</sub>)</th><th>log(SFR)<br>(<i>M</i><sub>☉</sub> yr<sup>−1</sup>)</th><th>log(sSFR)<br>(yr<sup>−1</sup>)</th><th>log(age)<br>(yr)</th><th><i>A</i><sub>V</sub><br>(mag)</th></tr></thead><tbody><tr><td>UVISTA-Y1</td><td>−22.48 ± 0.15</td><td>−1.5<sup>+0.4</sup><sub>−0.7</sub></td><td>0.45 ± 0.18</td><td>10.0<sup>+0.9</sup><sub>−0.4</sub></td><td>1.59<sup>+1.02</sup><sub>−9.55</sub></td><td>−8.4<sup>+1.8</sup><sub>−9.8</sub></td><td>7.30<sup>+1.42</sup><sub>−0.61</sub></td><td>0.9<sup>+0.0</sup><sub>−0.9</sub></td></tr><tr><td>UVISTA-Y2</td><td>−22.37 ± 0.20</td><td>−2.6<sup>+0.5</sup><sub>−0.5</sub></td><td>0.66 ± 0.21</td><td>9.0<sup>+0.3</sup><sub>−1.2</sub></td><td>1.98<sup>+0.65</sup><sub>−7.57</sub></td><td>−7.0<sup>+0.8</sup><sub>−8.6</sub></td><td>7.00<sup>+1.80</sup><sub>−0.50</sub></td><td>0.4<sup>+0.3</sup><sub>−0.4</sub></td></tr><tr><td>UVISTA-Y3a<sup>a</sup></td><td>−21.77 ± 0.32</td><td>−1.5<sup>+0.7</sup><sub>−0.9</sub></td><td>0.57 ± 0.39</td><td>9.8<sup>+1.3</sup><sub>−0.3</sub></td><td>−1.34<sup>+4.06</sup><sub>−7.01</sub></td><td>−11.1<sup>+4.9</sup><sub>−7.0</sub></td><td>8.00<sup>+0.80</sup><sub>−1.50</sub></td><td>0.0<sup>+1.1</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y3b<sup>a</sup></td><td>−21.23 ± 0.54</td><td>−3.2<sup>+2.2</sup><sub>−0.0</sub></td><td>−0.22 ± 1.82</td><td>8.7<sup>+0.1</sup><sub>−0.0</sub></td><td>−0.28<sup>+0.00</sup><sub>−0.07</sub></td><td>−9.0<sup>+0.0</sup><sub>−0.0</sub></td><td>7.40<sup>+0.00</sup><sub>−0.02</sub></td><td>0.0<sup>+0.0</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y3c<sup>a</sup></td><td>−21.37 ± 0.53</td><td>−2.0<sup>+1.7</sup><sub>−0.0</sub></td><td>0.81 ± 0.57</td><td>10.2<sup>+1.6</sup><sub>−0.5</sub></td><td>1.69<sup>+1.86</sup><sub>−28.23</sub></td><td>−8.5<sup>+2.3</sup><sub>−28.1</sub></td><td>8.70<sup>+0.10</sup><sub>−2.20</sub></td><td>0.8<sup>+1.1</sup><sub>−0.8</sub></td></tr><tr><td>UVISTA-Y4</td><td>−22.11 ± 0.24</td><td>−2.7<sup>+0.7</sup><sub>−0.4</sub></td><td>0.47 ± 0.26</td><td>9.9<sup>+0.5</sup><sub>−0.2</sub></td><td>1.23<sup>+0.61</sup><sub>−6.97</sub></td><td>−8.6<sup>+0.9</sup><sub>−7.0</sub></td><td>8.50<sup>+0.30</sup><sub>−3.1</sub></td><td>0.0<sup>+0.7</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y5</td><td>−22.34 ± 0.24</td><td>−1.7<sup>+0.8</sup><sub>−0.9</sub></td><td>0.63 ± 0.27</td><td>9.0<sup>+0.4</sup><sub>−1.1</sub></td><td>1.99<sup>+0.55</sup><sub>−7.68</sub></td><td>−7.0<sup>+0.8</sup><sub>−8.6</sub></td><td>7.30<sup>+1.40</sup><sub>−0.80</sub></td><td>0.4<sup>+0.2</sup><sub>−0.4</sub></td></tr><tr><td>UVISTA-Y6</td><td>−21.92 ± 0.26</td><td>−1.7<sup>+0.7</sup><sub>−0.8</sub></td><td>0.49 ± 0.32</td><td>9.7<sup>+1.1</sup><sub>−0.5</sub></td><td>1.36<sup>+1.33</sup><sub>−12.70</sub></td><td>−8.4<sup>+2.2</sup><sub>−12.9</sub></td><td>7.30<sup>+1.50</sup><sub>−0.80</sub></td><td>0.9<sup>+0.3</sup><sub>−0.9</sub></td></tr><tr><td>UVISTA-Y7</td><td>−21.74 ± 0.36</td><td>−2.0<sup>+0.7</sup><sub>−0.5</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>UVISTA-Y8</td><td>−21.76 ± 0.35</td><td>−2.8<sup>+0.9</sup><sub>−0.4</sub></td><td>0.91 ± 0.45</td><td>8.3<sup>+0.1</sup><sub>−1.4</sub></td><td>1.90<sup>+0.35</sup><sub>−1.41</sub></td><td>−6.4<sup>+0.2</sup><sub>−2.6</sub></td><td>6.50<sup>+2.29</sup><sub>−0.00</sub></td><td>0.0<sup>+0.5</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y9</td><td>−21.66 ± 0.34</td><td>−2.6<sup>+0.9</sup><sub>−0.6</sub></td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>UVISTA-Y10</td><td>−21.89 ± 0.31</td><td>−2.2<sup>+1.1</sup><sub>−0.7</sub></td><td>0.56 ± 0.39</td><td>8.3<sup>+0.0</sup><sub>−1.4</sub></td><td>1.80<sup>+0.47</sup><sub>−4.34</sub></td><td>−6.5<sup>+0.3</sup><sub>−5.6</sub></td><td>6.70<sup>+2.10</sup><sub>−0.20</sub></td><td>0.0<sup>+0.5</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y11</td><td>−22.04 ± 0.26</td><td>−1.8<sup>+0.5</sup><sub>−1.3</sub></td><td>0.67 ± 0.30</td><td>8.7<sup>+0.4</sup><sub>−1.2</sub></td><td>1.76<sup>+0.60</sup><sub>−7.62</sub></td><td>−7.0<sup>+0.8</sup><sub>−8.6</sub></td><td>7.30<sup>+1.44</sup><sub>−0.80</sub></td><td>0.3<sup>+0.2</sup><sub>−0.3</sub></td></tr><tr><td>UVISTA-Y12</td><td>−21.66 ± 0.40</td><td>−2.1<sup>+1.3</sup><sub>−0.8</sub></td><td>0.22 ± 0.64</td><td>9.1<sup>+0.9</sup><sub>−0.4</sub></td><td>0.17<sup>+2.22</sup><sub>−2.88</sub></td><td>−9.0<sup>+2.8</sup><sub>−3.1</sub></td><td>7.40<sup>+1.30</sup><sub>−0.90</sub></td><td>0.2<sup>+0.3</sup><sub>−0.2</sub></td></tr><tr><td>UVISTA-Y13</td><td>−21.39 ± 0.42</td><td>−1.0<sup>+0.7</sup><sub>−0.7</sub></td><td>0.50 ± 0.51</td><td>9.8<sup>+1.3</sup><sub>−0.3</sub></td><td>0.70<sup>+1.82</sup><sub>−9.08</sub></td><td>−9.1<sup>+1.8</sup><sub>−9.1</sub></td><td>7.50<sup>+1.28</sup><sub>−0.96</sub></td><td>0.8<sup>+0.3</sup><sub>−0.8</sub></td></tr><tr><td>UVISTA-Y14</td><td>−21.44 ± 0.40</td><td>−3.0<sup>+1.7</sup><sub>−0.1</sub></td><td>0.63 ± 0.54</td><td>9.3<sup>+1.2</sup><sub>−0.4</sub></td><td>0.52<sup>+2.03</sup><sub>−9.14</sub></td><td>−8.8<sup>+2.6</sup><sub>−9.3</sub></td><td>8.20<sup>+0.63</sup><sub>−1.70</sub></td><td>0.0<sup>+0.8</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y15</td><td>−21.50 ± 0.37</td><td>−2.7<sup>+0.9</sup><sub>−0.4</sub></td><td>−0.04 ± 0.68</td><td>8.8<sup>+0.2</sup><sub>−0.0</sub></td><td>−0.16<sup>+0.41</sup><sub>−0.02</sub></td><td>−9.0<sup>+0.6</sup><sub>−0.0</sub></td><td>7.40<sup>+0.01</sup><sub>−0.10</sub></td><td>0.0<sup>+0.0</sup><sub>−0.0</sub></td></tr><tr><td>UVISTA-Y16</td><td>−21.80 ± 0.33</td><td>−2.2<sup>+0.8</sup><sub>−0.8</sub></td><td>0.39 ± 0.40</td><td>8.6<sup>+0.3</sup><sub>−0.9</sub></td><td>1.62<sup>+0.56</sup><sub>−0.98</sub></td><td>−7.0<sup>+1.8</sup><sub>−1.8</sub></td><td>7.30<sup>+1.50</sup><sub>−0.80</sub></td><td>0.1<sup>+0.4</sup><sub>−0.1</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> These 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 \pm 0.16$ mag, $\beta = -1.8 \pm 0.7$, $u-g = 0.58 \pm 0.16$ mag, $\log(M_*/M_{\odot}) = 9.9_{-0.3}^{+0.6}$, $\log(SFR/M_{\odot}/\text{yr}^{-1}) = 1.63_{-3.77}^{+0.38}$, $\log(\text{sSFR}/\text{yr}^{-1}) = -8.2_{-3.8}^{+0.9}$, $\log(\text{age}/\text{yr}) = 8.20_{-1.16}^{+0.60}$, and $A_V = 0.5_{-0.5}^{+0.5}$ mag.

<sup>b</sup> After 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.

<table><caption>**Table 5**
Observed and Rest-frame Properties for Candidate *z* ~ 8 Galaxies Identified in the UltraVISTA DR3 Observations</caption><thead><tr><th>Quantity</th><th>25%</th><th>Median</th><th>75%</th><th>25% Uncertainties</th><th>Median Uncertainties</th><th>75% Uncertainties</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td>8.05</td><td>8.40</td><td>8.62</td><td>+0.60/-0.61</td><td>+0.69/-0.73</td><td>+0.96/-1.15</td></tr><tr><td>M<sub>UV</sub> [mag]</td><td>-22.0</td><td>-21.8</td><td>-21.6</td><td>±0.3</td><td>±0.3</td><td>±0.4</td></tr><tr><td>UV β</td><td>-2.68</td><td>-2.17</td><td>-1.73</td><td>+0.70/-0.40</td><td>+0.77/-0.65</td><td>+1.01/-0.79</td></tr><tr><td>(u-g)<sub>rest</sub> [mag]</td><td>0.42</td><td>0.53</td><td>0.65</td><td>±0.29</td><td>±0.39</td><td>±0.55</td></tr><tr><td>log(M<sub>*</sub>/M<sub>☉</sub>)</td><td>8.71</td><td>9.07</td><td>9.76</td><td>+0.32/-0.24</td><td>+0.46/-0.44</td><td>+1.14/-1.15</td></tr><tr><td>M<sub>*</sub>/L<sub>UV</sub> [M<sub>☉</sub>/L<sub>☉</sub>]</td><td>0.005</td><td>0.010</td><td>0.044</td><td>+0.008/-0.004</td><td>+0.015/-0.014</td><td>+0.034/-0.098</td></tr><tr><td>M<sub>*</sub>/L<sub>U</sub> [M<sub>☉</sub>/L<sub>☉</sub>]</td><td>0.013</td><td>0.048</td><td>0.101</td><td>+0.026/-0.010</td><td>+0.038/-0.056</td><td>+0.098/-0.244</td></tr><tr><td>M<sub>*</sub>/L<sub>g</sub> [M<sub>☉</sub>/L<sub>☉</sub>]</td><td>0.017</td><td>0.064</td><td>0.133</td><td>+0.026/-0.007</td><td>+0.043/-0.089</td><td>+0.078/-0.188</td></tr><tr><td>log(SFR/M<sub>*</sub>/\text{yr}^{-1})</td><td>0.3</td><td>1.5</td><td>1.8</td><td>+0.5/-2.1</td><td>+0.6/-7.3</td><td>+1.8/-9.1</td></tr><tr><td>log(sSFR/\text{yr}^{-1})</td><td>-9.0</td><td>-8.4</td><td>-7.0</td><td>+0.7/-2.9</td><td>+0.9/-7.8</td><td>+2.5/-9.2</td></tr><tr><td>log(age/\text{yr})</td><td>7.30</td><td>7.35</td><td>7.75</td><td>+0.47/-0.35</td><td>+1.35/-0.80</td><td>+1.50/-1.13</td></tr><tr><td>A<sub>V</sub> [mag]</td><td>0.00</td><td>0.15</td><td>0.60</td><td>+0.22/-0.00</td><td>+0.32/-0.15</td><td>+0.56/-0.60</td></tr></tbody></table>

**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 $\beta$ were measured following the procedure described in Section 6.2. The last two columns present the first and third quartiles of uncertainties, respectively.

(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).

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$\beta$ 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
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The best-fit stellar masses with emission lines included compared to those derived with the standard BC03 models without emission lines. The latter masses (where the models ignore line emission) are higher by ~0.43 dex on average, consistent with the results of Labbé et al. (2013), with individual galaxies differing by up to 1 dex. One might expect more accurate masses from standard BC03 models if one excludes the 4.5 µm band (contaminated by [O III]+Hβ emission) when performing the fitting, but the estimated stellar masses are still found to be 0.23 dex higher on average. This mismatch between the BC03 model fit results (without the emission lines) and the fit results with emission lines included may be due to the contribution of the [O II] line to the 3.6 µm band flux measurements. From the present exercise, we can see how important it is to fully consider nebular emission when estimating stellar population parameters.</figcaption></figure>

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 \sim 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]$_{5007}$ 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]$_{3727}$) 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 \gtrsim 8$ galaxies.

The typical estimated stellar masses for bright sources in our $z \sim 8$ selection (see Table 5) are $10^{9.1^{+0.5}_{-0.4}} M_{\odot}$, with the SFRs of $32^{+44}_{-32} M_{\odot} \text{yr}^{-1}$, specific SFR of $4^{+8}_{-4} \text{Gyr}^{-1}$, stellar ages of $\sim 22^{+69}_{-22} \text{Myr}$, and low dust content $A_V = 0.15^{+0.30}_{-0.15} \text{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$_{1600}$ 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.01 M_{\odot}/L_{\odot}$.

We find a positive although marginal correlation of the $M_{\star}/L_{UV,1600}$ with the UV slope for our $z \sim 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_{\star}/L_{\odot}$ ratios, such that redder galaxies exhibit higher $M_{\star}/L_{\odot}$, 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 \sim 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 \sim 8$ sample do not correspond to higher $M_{\star}/L_{\odot}$. Instead, the optically reddest galaxies tend to have the lowest $M_{\star}/L_{\odot}$. 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 $\beta$ and $u-g$ in Figure 10.

### 6.5. Volume Density of Bright $z \sim 8$ and $z \sim 9$ Galaxies

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 \sim 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} \sim 9$, we also explored the contribution to the UV LF at $z \sim 9$ from the five sources with $z_{phot} \ge 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 \sim 9$ LF, bringing to six the total number of sources used for the $z \sim 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 \sim 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 $\propto(1+z)^{-1}$ (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 $\beta = -1.8 \pm 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_{\nu}$ and strong emission lines with fixed
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Left panel: the rest-frame UV-continuum slope $\beta$ vs. the best-fit stellar mass-to-UV light ratio. The dashed line marks a tentative linear correlation. Right panel: the rest-frame $u-g$ color vs. the best-fit stellar mass-to-optical light ratio. At low to intermediate redshift $z \sim 2$, a tight relation exists between rest-frame $u-g$ colors and $M_{*}/L_{g}$ ratios with unity slope, such that redder galaxies exhibit higher $M_{*}/L_{g}$. The gray dashed line shows the relation derived by Szomoru et al. (2013) at $z \sim 2$. The orange dashed curve shows the relation for our BC03 models, including emission lines. The orange triangles mark the age of the stellar population, starting from log(age/yr) = 6 to log(age/yr) = 8.5, in steps of 0.5 dex. While we find a marginal positive correlation with UV slope for our $z \sim 8$ sample, there is no clear relation between $u-g$ colors and $M_{*}/L_{g}$. Instead, the optically reddest galaxies tend to have the lowest $M_{*}/L_{g}$ ratios. This likely reflects the impact of strong emission lines on the $g$-band fluxes.</figcaption></figure>

rest-frame EW(H$\alpha$+[N II]+[S II]) = 300 Å and rest-frame EW([O III]+H$\beta$) = 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 \sim 24.0-24.5$ mag and $24.5-25.0$ is $5.3 \times 10^6$ Mpc$^3$ and $2.6 \times 10^6$ Mpc$^3$, respectively.

We estimate constraints on the bright end of the UV LF, adopting the $V_{\text{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 \sim 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_{\text{UV}} \le -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 \ge 8.6$ (after forcing the inclusion of UVISTA-Y5), we considered these galaxies to belong to the $z \sim 9$ redshift bin and computed the associated number densities accordingly.

The resulting $z \sim 9$ LF is presented in the bottom panel of Figure 14 and in Table 6.

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 \sim 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 \sim 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 \sim 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 \sim 8$ candidates have apparent magnitudes fainter than $H \sim 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 \sim 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 \sim 8$ LF is consistent with the previous measurement at $M_{\text{UV}} \le -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 ($\sim L^*$) LBGs. However, the revised analysis performed for the current work showed that one of the sources
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Top panel: The blue points with error bars mark our estimates of volume density associated to the sample of candidate luminous $z \sim 8$ galaxies considered in this work. For comparison, we also present recent UV LF determinations at $z \sim 8$ from empty field studies, as indicated by the legend (we arbitrarily shifted the measurement of Bridge et al. 2019 by +0.05 mag to improve readability). Bottom panel: Here we compare our $z \sim 9$ volume density estimate from the five sources with $z_{phot} \ge 8.6$ (blue points) to measurements of the UV LF at $z \sim 9$. The blue open circle at the bright end corresponds to the measurement obtained forcing UVISTA-Y5 into the $z \sim 9$ sample (see Section 6.5 for details), while the blue open circle at $M_{UV} = -21.2$ mag marks the volume density measurement for the faintest luminosity bin, where our sample is likely incomplete. The magenta curve presents the bright end of the dual power law from Bowler et al. (2017) evolved to $z \sim 9$ following Bouwens et al. (2016), and whose characteristic density has been adjusted to match that of the Schechter function at the characteristic luminosity.</figcaption></figure>

previously considered to be at $z_{phot} \sim 7.5$ is actually at $z_{phot} \gtrsim 8$, thus increasing its luminosity and balancing the final volume density.

For $M_{UV} \lesssim -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\sigma$.

At $M_{UV} \gtrsim -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)

<table><caption>**Table 6**
$V_{max}$ Determinations of the UV LF</caption><thead><tr><th>M<sub>UV</sub><br>(mag)</th><th>φ<br>(×10<sup>-6</sup> mag<sup>-1</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td colspan="2">z ~ 8</td></tr><tr><td>-22.55</td><td>0.76<sup>+0.74</sup><sub>-0.41</sub></td></tr><tr><td>-22.05</td><td>1.38<sup>+1.09</sup><sub>-0.66</sub></td></tr><tr><td>-21.55<sup>a</sup></td><td>4.87<sup>+2.01</sup><sub>-1.41</sub></td></tr><tr><td colspan="2">z ~ 9</td></tr><tr><td>-22.35<sup>b</sup></td><td>0.43<sup>+0.99</sup><sub>-0.36</sub></td></tr><tr><td>-22.00</td><td>0.43<sup>+0.98</sup><sub>-0.36</sub></td></tr><tr><td>-21.60</td><td>1.14<sup>+1.50</sup><sub>-0.73</sub></td></tr><tr><td>-21.20<sup>a</sup></td><td>1.64<sup>+2.16</sup><sub>-1.06</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> This 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.

<sup>b</sup> The volume density in this luminosity bin was obtained forcing UVISTA-Y5 into the sample of galaxies at $z_{phot} \ge 8.6$ (i.e., our $z \sim 9$ LBG sample). Its nominal $z_{phot} = 8.596$ would exclude it from the sample of $z \sim 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.

presented the $z \sim 8$ LF from eight $M_{UV} \gtrsim -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 $\sim 5 \times$ higher than what we estimate for our sample, and their measurements are only consistent at $\sim 3\sigma$. 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} \sim 8.4$).

In the lower panel of Figure 14 we present our estimates of the $z \sim 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} \lesssim -22$ mag our new bright $z \sim 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} \sim -22$ mag is consistent at $1\sigma$ 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 $\sim 2\sigma$, our density being $\sim 9 \times$ 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 \sim 9$ sources have available *Spitzer/IRAC* or deep Y098 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 \sim 9$, applying the relations of Bouwens et al. (2016) to the DPL found at $z \sim 7$ by Bowler et al. (2017,
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Step-wise LF at $z \sim 8$ obtained combining the $V_{max}$ estimates at the bright end combined to those of Bouwens et al. (2015) at $M_{UV} > -22$ mag. The discrete LF measurements can be represented by either a Schechter (1976, blue solid curve) or double-power-law (red solid curve) form, with a marginal preference for this latter model. We also show the halo-mass function (HMF) at $z \sim 8$ (orange dashed line) scaled by a fixed $M_{halo}/L_{UV}$ to match the knee of our derived UV LF at $z \sim 8$. The high-mass-end slope of the HMF is similar to the effective slope of the UV LF at the bright end. The difference between the low-mass-end slope of the HMF and the faint-end slope of the LF is $\Delta\alpha \sim 0.3$.</figcaption></figure>

see also Stefanon et al. 2017b). Indeed, the excess in number density we observe for $M_{UV} < -22$ mag at $z \sim 9$ (introduced by forcing UVISTA-Y5 into the $z \sim 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 \sim 9$, which instead needs larger samples. We will discuss the shape of the $z \sim 8$ LF in more detail in Section 6.7.

## 6.6. Combination of Present Constraints with Faint $z \sim 8$ LF Results

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 $\alpha$ and turn over magnitude $M^*$, both bright and faint objects are needed to constrain $\alpha$ and $M^*$. The similar redshift distributions expected for bright galaxies selected by our $z \sim 8$ criteria and those selected in the fainter Bouwens et al. (2015) samples (see Figure 2) make it possible to combine our $z \sim 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 \sim 8$ is presented in Figure 15. We determined the Schechter function parameters $M^*$, $\alpha$, and $\phi^*$ minimizing the $\chi^2$, and obtaining $\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.

Our sample of bright $z \sim 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 $\sim 1.2\sigma$, while the faint-end slope $\alpha$ is consistent at $1\sigma$.

In Figure 16 we also compare our estimated Schechter parameters to their evolution over a wide range of redshift,

$4 \lesssim z \lesssim 10$ from Bouwens et al. (2015). Our result confirms the picture of marginal evolution of $M^*$ for $z \gtrsim 3-4$, but significant evolution of $\alpha$ and $\phi^*$. This conclusion was first drawn by Bouwens et al. (2015) using LF results from $z \sim 7$ to $z \sim 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 \sim 0.2$ from $z \sim 7$ to $z \sim 5$.

## 6.7. The Shape of the LF at $z \sim 8$

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 \sim 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 \sim 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):

$$ \phi(M) = \frac{\phi^*}{10^{0.4(\alpha+1)(M-M^*)} + 10^{0.4(\beta+1)(M-M^*)}} $$

where $\alpha$ and $\beta$ are the faint-end and bright-end slopes, respectively, $M^*$ is the transition luminosity between the two power-law regimes, and $\phi^*$ is the normalization.

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

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 \lesssim 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 \sim 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,
<figure><img src="image_19.png" /><figcaption>**Figure 16.** 68% and 95% likelihood contours on the Schechter parameters at $z \sim 8$ (red) derived in the present work shown relative to the Schechter parameter estimates found for the LF results at $z \sim 4$-to-10 by Bouwens et al. (2015), as specified by the legend. The observations seem to point toward a clear increase in $\phi^*$ and flattening of $\alpha$ with cosmic time.</figcaption></figure>

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 \sim 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 \sim 8$ constraints, we find $M^* = -20.81$ mag, $\sim 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} \lesssim -22.5$ mag) of the $z \sim 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<fnref n="16" /> (Murray et al. 2013), assuming a Planck Collaboration et al. (2016) cosmology, with $\Omega_m = 0.2678$, $\Omega_b = 0.049$, $H_0 = 67.04$, $n_s = 0.962$, and $\sigma_8 = 0.8347$.

The scaled HMF looks similar to the UV LF at $z \sim 8$. We observe only a small $\Delta\alpha \sim 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 ($\Delta\beta \sim 0.2–0.3$). Interestingly, the bright end of the double power law overlaps with the HMF for $M_{UV} \lesssim -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 \sim 9$, due to the challenges in determining the total number of bright $z \sim 9$ galaxies over UltraVISTA. Clearly, if any significant number of the candidates do prove to be bona fide $z \sim 9$ galaxies, they would favor more of a power-law form to the bright end of the $z \sim 9$ LF.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** The evolution of the UV luminosity density of galaxies at $4 < z < 9$ (filled blue circle) brightward of $-21.5$ mag using the present search over COSMOS/UltraVISTA. The green squares mark the luminosity density of galaxies obtained from the Schechter parameterization by Bouwens et al. (2015) at $z \sim 4–10$. Also shown are the recent estimates at $z \sim 9$ of Bouwens et al. (2016) and at $z \sim 10$ of Oesch et al. (2018). The shaded region shows a linear fit to the evolutionary trend to $z \sim 10$ in the UV luminosity density preferred at 68% confidence.</figcaption></figure>

## 6.8. Evolution of the UV Luminosity Density for Luminous Galaxies from $z \sim 10$ to $z \sim 4$

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., $\sim -21.5$ mag) at $z \sim 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 \sim 8$ is $10^{23.87_{-0.68}^{+0.58}}$ ergs/s/Hz/Mpc³ (68% c.l.).

In Figure 17 we also present the luminosity density determinations and uncertainties that we derive integrating to

<footnote n="16">¹⁶ http://hmf.icrar.org/hmf finding form/create/</footnote>
$M_{\text{UV}} = -21.5$ mag the recent LF results of Bouwens et al. (2015) at $z \sim 4$, $z \sim 5$, $z \sim 6$, $z \sim 7$, and $z \sim 8$; the results of Bouwens et al. (2016) at $z \sim 9$; and the results at $z \sim 10$ from Oesch et al. (2018), together with our best-fit constraints and $1\sigma$ 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 \sim 8$ LF of Bouwens et al. (2015) and with the upper limit at $z \sim 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 \sim 10$ are $\sim 1.7\sigma$ and $\sim 1.1\sigma$ 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 $\sim 500$ Myr of cosmic history (see e.g., the extensive discussion in Oesch et al. 2018).

## 7. Conclusions

Using deep infrared data from the COSMOS/UltraVISTA program, we have identified 16 new ultrabright $H \sim 24.8-25.6$ mag galaxy candidates at $z \sim 8$. The new candidates are among the brightest yet found at these redshifts, $\gtrsim 0.5$ mag brighter than found over CANDELS, providing improved constraints at the bright end of the UV LF, and providing excellent targets for follow-up at longer wavelengths and with spectroscopy.

The spectral slope of the UV-continuum $\beta$, parameterized as $f_{\lambda} \propto \lambda^{\beta}$; (Meurer et al. 1999) for the bright $z \sim 8$, sample is $\beta = -2.2 \pm 0.6$, which is bluer but still consistent with the UV-continuum slopes found for bright $-22 < M_{\text{UV}} < -21$ galaxies at $z = 6$ ($\beta = -1.55 \pm 0.17$) and $z \sim 7$ ($\beta = -1.75 \pm 0.18$) by Bouwens et al. (2014) and suggests a continuing trend toward bluer $\beta$s at higher redshifts. The typical estimated stellar masses for bright sources in our $z \sim 8$ selection are $10^{9.1^{+0.5}_{-0.4}} M_{\odot}$, with the SFRs of $32_{-32}^{+44} M_{\odot} \text{yr}^{-1}$, specific SFR of $4_{-4}^{+8} \text{Gyr}^{-1}$, stellar ages of $\sim 22_{-22}^{+69}$ Myr, and low dust content $A_{\text{V}} = 0.15_{-0.15}^{+0.30}$ mag, with the properties of individual galaxies spanning a large range of values.

Using public catalogs we checked the lensing magnification from close, lower-redshift sources. We find that four sources are likely subject to magnifications of approximately $1.5 \times$. Nevertheless, the effect on the UV LF is marginal.

We use the candidate galaxies to constrain the bright end of the $z \sim 8$ UV LF. Combining our ultrabright sample with candidates found over CANDELS, HUDF, and HFF field data allows us to constrain on the $z \sim 8$ LF. Assuming a Schechter function, the best-fit characteristic magnitude is $M * (z = 8) = -20.95_{-0.35}^{+0.30}$ mag with a very steep faint-end slope $\alpha = -2.15_{-0.19}^{+0.20}$. Our $z \sim 8$ LF results can be equally well represented adopting a functional form where the effective slope is steeper at the bright end of the LF than at the faint end, such as for a DPL. Our results rule out the use of a single power-law in representing the $z \sim 8$ LF.

We note that, despite much recent progress, the lack of spectra and deep high-resolution imaging still limit us in establishing the reliability of high-redshift galaxy selections, in particular for rare luminous galaxies that constrain the bright end of the UV luminosity and mass functions where any contamination has a very large impact. While care is taken in estimating the completeness and contamination rates, these still rely on assumed SEDs. Ultimately, spectroscopy is needed

to validate these assumptions. While recent results suggest ALMA as a potentially efficient machine for the study of emission lines (e.g., Smit et al. 2018), currently this is still hard in the rest-frame UV and optical, due to long integration times, low multiplexing, and the reduced observable emission of Ly$\alpha$ likely caused by the increasing neutral hydrogen fraction $z > 6$ (e.g., Schenker et al. 2014), but will be possible in the future with JWST and the next generation of extremely large ground-based telescopes.

The authors would like to thank the referee for a careful reading and for the many helpful, positive, and very constructive comments that helped improving the quality of the paper. The authors wish to acknowledge the UltraVISTA team who conducted and followed the observations and image processing. The authors are thankful to the COSMOS collaboration for their continued efforts in making available part of the data used in this work. KIC acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. The Cosmic Dawn Center is funded by the Danish National Research Foundation. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan. The HSC collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University. This research has benefitted from the SpeX Prism Spectral Libraries, maintained by Adam Burgasser at http://pono.ucsd.edu/~adam/browndwarfs/spexprism.

# Appendix A

## Monte Carlo Assessment of the Multi-component Nature of UVISTA-Y3

Because the three components of UVISTA-Y3 are characterized by low S/N on the DASH footprint adopted for their identification (S/N ~ 4.5, 2.9 and ~2.2), one might wonder whether the detected splitting into three components is instead the result of background noise acting on a single, extended source.

To test this hypothesis we implemented the following Monte Carlo procedure. We generated a table of 20 random positions
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Each panel presents one of the random realizations of the intrinsic (i.e., before being convolved with the WFC3 H<sub>160</sub>-band PSF and added to the DASH mosaic) exponential disk created to test the multi-component nature of UVISTA-Y3. The side of each stamp is 2"/0.</figcaption></figure>

on the footprint of the DASH mosaic with similar background noise properties. For each one of these positions we created an elongated disk, with minor-to-major axis ratio *b*/*a* drawn from a pool of random values 0.05 < *b*/*a* < 0.25 with an exponential luminosity profile with effective radius *r*<sub>e</sub> = 1 kpc, consistent with recent rest-frame UV size estimates for bright LBGs at *z* > 6 (e.g., Oesch et al. 2016; Bowler et al. 2017; Stefanon et al. 2017b). The choice of low values for *b*/*a* was guided by the relatively large separation between the three components, which is difficult to obtain when more compact morphologies are considered. The total flux density of each exponential disk was set to be equal to the *H*-band flux of UVISTA-Y3 when considered as a single source. Each exponential disk was then convolved with the WFC3 H<sub>160</sub>-band PSF, randomly rotated and added to the original DASH image.

In Figure 18 we present the 20 random realizations of the exponential disk, before convolution with the WFC3 PSF, while in Figure 19 we present the image stamps of the DASH mosaics after the synthetic exponential disks have been added.

Based on simple visual inspection, we see no indication in these simulated images for a multiple component structure. Finally we run SExtractor using the same set of parameters adopted for the original deblending and found that none of the synthetic sources were split into two or more components. This test therefore increased our confidence in the multiple nature of UVISTA-Y3. It is worth remarking that as a result of the low S/N of the deblended photometry for each component, there are substantially larger uncertainties in the derived physical parameters for each component.

We complemented this first assessment with a second Monte Carlo simulation in which we adopt an effective radius $r_e = 3$ kpc, similar to the sizes of luminous high-redshift LBGs when potentially multiple sources are considered as a single object (e.g., Bowler et al. 2017), and $0.1 < b/a < 0.7$. We present the result of this simulation in Figures 20 and 21. Using a procedure similar to that applied for the $r_e = 1$ kpc case, we do not find evidence for multiple components even when large $r_e$ are considered, increasing the confidence in our interpretation of the three components in UVISTA-Y3.
<figure><img src="image_22.png" /><figcaption>**Figure 19.** Each image stamp (2''0 side) is centered on the corresponding synthetic exponential disk of Figure 18, after being convolved with the WFC3 PSF and added to the DASH mosaic at locations with noise properties similar to those where UVISTA-Y3 lies. In the central panel of the bottom row, the object on the right is a real source present in the observations (and therefore bears no relation to the Monte Carlo simulations we perform). None of the simulated sources show apparent multi-component structure (as UVISTA-Y3 seems to show) as a result of noise in the background.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20. Similar to Figure 18 but with $r_c = 3$ kpc.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 21.** Similar to Figure 19 but with $r_c = 3$ kpc. The image stamps have been smoothed with a Gaussian of 0.1 FWHM (equal to the value used for Figure 4) to improve the contrast.</figcaption></figure>
## Appendix B
Flux Density Estimates

In Tables 7–9 we list the flux density estimates and associated 1σ uncertainties for the full sample of $z \sim 8$ candidate galaxies presented in this work.

<table><caption>Table 7
Flux Density Measurements for UVISTA-Y1 to UVISTA-Y5</caption><thead><tr><th>Filter<br>Name</th><th>UVISTA-Y1<br>(nJy)</th><th>UVISTA-Y2<br>(nJy)</th><th>UVISTA-Y3a<br>(nJy)</th><th>UVISTA-Y3b<br>(nJy)</th><th>UVISTA-Y3c<br>(nJy)</th><th>UVISTA-Y4<br>(nJy)</th><th>UVISTA-Y5<br>(nJy)</th></tr></thead><tbody><tr><td>CFHTLS <i>u</i>*</td><td>...</td><td>-6 ± 15</td><td>7 ± 15</td><td>-7 ± 15</td><td>-5 ± 14</td><td>-10 ± 15</td><td>-13 ± 17</td></tr><tr><td>SSC <i>B</i></td><td>-2 ± 7</td><td>-1 ± 10</td><td>1 ± 7</td><td>-6 ± 8</td><td>-10 ± 9</td><td>5 ± 9</td><td>4 ± 9</td></tr><tr><td>HSC <i>g</i></td><td>8 ± 15</td><td>2 ± 17</td><td>-22 ± 16</td><td>5 ± 17</td><td>2 ± 17</td><td>-8 ± 17</td><td>-10 ± 20</td></tr><tr><td>CFHTTLS <i>g</i></td><td>...</td><td>-1 ± 13</td><td>-10 ± 12</td><td>-11 ± 12</td><td>-6 ± 12</td><td>2 ± 13</td><td>3 ± 15</td></tr><tr><td>SSC <i>V</i></td><td>-12 ± 17</td><td>-14 ± 22</td><td>-5 ± 21</td><td>0 ± 21</td><td>-7 ± 20</td><td>29 ± 20</td><td>-1 ± 24</td></tr><tr><td>HSC <i>r</i></td><td>-7 ± 14</td><td>-4 ± 17</td><td>-5 ± 14</td><td>-5 ± 15</td><td>1 ± 16</td><td>7 ± 17</td><td>-3 ± 18</td></tr><tr><td>CFHTTLS <i>r</i></td><td>...</td><td>-12 ± 19</td><td>-1 ± 17</td><td>-14 ± 17</td><td>14 ± 19</td><td>7 ± 20</td><td>-15 ± 23</td></tr><tr><td>SSC <i>r</i><sup>+</sup></td><td>7 ± 16</td><td>-15 ± 20</td><td>1 ± 17</td><td>5 ± 17</td><td>10 ± 18</td><td>10 ± 18</td><td>-29 ± 23</td></tr><tr><td>SSC <i>i</i><sup>+</sup></td><td>-9 ± 22</td><td>-139 ± 89</td><td>-17 ± 27</td><td>-27 ± 27</td><td>-21 ± 28</td><td>-103 ± 90</td><td>-36 ± 26</td></tr><tr><td>CFHTTLS <i>y</i></td><td>...</td><td>-18 ± 25</td><td>-13 ± 22</td><td>17 ± 20</td><td>2 ± 24</td><td>-14 ± 25</td><td>-11 ± 28</td></tr><tr><td>CFHTTLS <i>i</i></td><td>...</td><td>-28 ± 26</td><td>-20 ± 22</td><td>-19 ± 21</td><td>-11 ± 24</td><td>3 ± 26</td><td>-11 ± 29</td></tr><tr><td>HSC <i>i</i></td><td>21 ± 21</td><td>-44 ± 24</td><td>-12 ± 21</td><td>-21 ± 20</td><td>-2 ± 24</td><td>-10 ± 24</td><td>-20 ± 27</td></tr><tr><td>CFHTTLS <i>z</i></td><td>...</td><td>-13 ± 55</td><td>-22 ± 46</td><td>0 ± 45</td><td>-15 ± 52</td><td>7 ± 54</td><td>6 ± 63</td></tr><tr><td>HSC <i>z</i></td><td>9 ± 31</td><td>-1 ± 35</td><td>1 ± 30</td><td>-30 ± 29</td><td>-51 ± 35</td><td>-50 ± 35</td><td>-27 ± 39</td></tr><tr><td>SSC <i>z</i><sup>+</sup></td><td>51 ± 64</td><td>-102 ± 75</td><td>19 ± 72</td><td>-52 ± 75</td><td>79 ± 72</td><td>43 ± 71</td><td>-51 ± 93</td></tr><tr><td>HSC <i>y</i></td><td>-31 ± 76</td><td>31 ± 85</td><td>-5 ± 81</td><td>-86 ± 80</td><td>-83 ± 84</td><td>33 ± 84</td><td>-91 ± 98</td></tr><tr><td>UVISTA <i>Y</i></td><td>18 ± 48</td><td>29 ± 53</td><td>55 ± 50</td><td>-12 ± 49</td><td>80 ± 53</td><td>173 ± 54</td><td>-42 ± 68</td></tr><tr><td>UVISTA <i>J</i></td><td>324 ± 50</td><td>410 ± 61</td><td>120 ± 57</td><td>111 ± 57</td><td>33 ± 59</td><td>432 ± 65</td><td>235 ± 66</td></tr><tr><td>UVISTA <i>H</i></td><td>455 ± 61</td><td>432 ± 78</td><td>229 ± 68</td><td>135 ± 67</td><td>145 ± 70</td><td>392 ± 86</td><td>393 ± 86</td></tr><tr><td>UVISTA <i>K</i><sub>S</sub></td><td>480 ± 77</td><td>275 ± 86</td><td>197 ± 77</td><td>68 ± 72</td><td>110 ± 82</td><td>266 ± 110</td><td>321 ± 102</td></tr><tr><td>IRAC 3.6 µm</td><td>623 ± 85</td><td>492 ± 50</td><td>220 ± 78</td><td>50 ± 80</td><td>214 ± 79</td><td>620 ± 68</td><td>289 ± 74</td></tr><tr><td>IRAC 4.5 µm</td><td>931 ± 109</td><td>799 ± 57</td><td>389 ± 99</td><td>-15 ± 101</td><td>393 ± 90</td><td>682 ± 108</td><td>589 ± 86</td></tr><tr><td>IRAC 5.8 µm</td><td>-2893 ± 2568</td><td>688 ± 1702</td><td>-1028 ± 3081</td><td>4174 ± 3169</td><td>2700 ± 2979</td><td>-1686 ± 1819</td><td>-1978 ± 4831</td></tr><tr><td>IRAC 8.0 µm</td><td>1423 ± 3021</td><td>1384 ± 2105</td><td>-3418 ± 4231</td><td>-753 ± 4341</td><td>-1036 ± 3907</td><td>-795 ± 2123</td><td>499 ± 6310</td></tr></tbody></table>

<table><caption>Table 8
Flux Density Measurements for UVISTA-Y6 to UVISTA-Y12</caption><thead><tr><th>Filter<br>Name</th><th>UVISTA-Y6<br>(nJy)</th><th>UVISTA-Y7<br>(nJy)</th><th>UVISTA-Y8<br>(nJy)</th><th>UVISTA-Y9<br>(nJy)</th><th>UVISTA-Y10<br>(nJy)</th><th>UVISTA-Y11<br>(nJy)</th><th>UVISTA-Y12<br>(nJy)</th></tr></thead><tbody><tr><td>CFHTTLS <i>u</i>*</td><td>8 ± 14</td><td>-10 ± 15</td><td>-15 ± 15</td><td>...</td><td>-10 ± 14</td><td>-3 ± 16</td><td>...</td></tr><tr><td>SSC <i>B</i></td><td>-11 ± 9</td><td>-12 ± 9</td><td>0 ± 8</td><td>4 ± 10</td><td>0 ± 9</td><td>3 ± 7</td><td>-12 ± 9</td></tr><tr><td>HSC <i>g</i></td><td>4 ± 16</td><td>-13 ± 19</td><td>0 ± 17</td><td>-2 ± 17</td><td>-3 ± 17</td><td>-19 ± 15</td><td>1 ± 16</td></tr><tr><td>CFHTTLS <i>g</i></td><td>1 ± 13</td><td>-17 ± 13</td><td>-5 ± 13</td><td>...</td><td>-21 ± 13</td><td>6 ± 12</td><td>...</td></tr><tr><td>SSC <i>V</i></td><td>0 ± 21</td><td>11 ± 21</td><td>2 ± 19</td><td>-1 ± 20</td><td>11 ± 19</td><td>-1 ± 17</td><td>-1 ± 20</td></tr><tr><td>HSC <i>r</i></td><td>11 ± 15</td><td>2 ± 21</td><td>16 ± 17</td><td>3 ± 16</td><td>9 ± 16</td><td>-3 ± 14</td><td>-13 ± 16</td></tr><tr><td>CFHTTLS <i>r</i></td><td>8 ± 21</td><td>-13 ± 23</td><td>-5 ± 19</td><td>...</td><td>5 ± 19</td><td>1 ± 18</td><td>...</td></tr><tr><td>SSC <i>r</i><sup>+</sup></td><td>20 ± 19</td><td>-34 ± 20</td><td>-21 ± 18</td><td>-8 ± 19</td><td>-14 ± 18</td><td>-1 ± 17</td><td>-18 ± 20</td></tr><tr><td>SSC <i>i</i><sup>+</sup></td><td>2 ± 21</td><td>-42 ± 32</td><td>-26 ± 29</td><td>-24 ± 30</td><td>-31 ± 29</td><td>-29 ± 27</td><td>-2 ± 31</td></tr><tr><td>CFHTTLS <i>y</i></td><td>22 ± 26</td><td>-14 ± 25</td><td>-5 ± 25</td><td>...</td><td>-5 ± 25</td><td>14 ± 22</td><td>...</td></tr><tr><td>CFHTTLS <i>i</i></td><td>5 ± 29</td><td>-16 ± 27</td><td>-11 ± 26</td><td>...</td><td>12 ± 25</td><td>8 ± 22</td><td>...</td></tr><tr><td>HSC <i>i</i></td><td>1 ± 23</td><td>8 ± 30</td><td>-7 ± 24</td><td>-3 ± 24</td><td>-18 ± 24</td><td>14 ± 21</td><td>-16 ± 24</td></tr><tr><td>CFHTTLS <i>z</i></td><td>-13 ± 59</td><td>-18 ± 55</td><td>-5 ± 55</td><td>...</td><td>27 ± 54</td><td>-13 ± 46</td><td>...</td></tr><tr><td>HSC <i>z</i></td><td>17 ± 33</td><td>10 ± 40</td><td>-16 ± 35</td><td>-21 ± 35</td><td>50 ± 34</td><td>43 ± 32</td><td>14 ± 34</td></tr><tr><td>SSC <i>z</i><sup>+</sup></td><td>69 ± 85</td><td>19 ± 77</td><td>16 ± 73</td><td>-74 ± 78</td><td>54 ± 72</td><td>-75 ± 63</td><td>-41 ± 76</td></tr><tr><td>HSC <i>y</i></td><td>89 ± 80</td><td>62 ± 88</td><td>55 ± 85</td><td>33 ± 87</td><td>-8 ± 85</td><td>91 ± 83</td><td>14 ± 84</td></tr><tr><td>UVISTA <i>Y</i></td><td>16 ± 51</td><td>-9 ± 57</td><td>-31 ± 54</td><td>101 ± 53</td><td>20 ± 53</td><td>-17 ± 45</td><td>-79 ± 53</td></tr><tr><td>UVISTA <i>J</i></td><td>211 ± 53</td><td>189 ± 64</td><td>192 ± 61</td><td>204 ± 59</td><td>242 ± 61</td><td>158 ± 55</td><td>164 ± 58</td></tr><tr><td>UVISTA <i>H</i></td><td>280 ± 66</td><td>231 ± 76</td><td>241 ± 77</td><td>246 ± 78</td><td>275 ± 79</td><td>297 ± 71</td><td>207 ± 76</td></tr><tr><td>UVISTA <i>K</i><sub>S</sub></td><td>271 ± 82</td><td>204 ± 84</td><td>85 ± 87</td><td>116 ± 87</td><td>211 ± 88</td><td>201 ± 77</td><td>191 ± 83</td></tr><tr><td>IRAC 3.6 µm</td><td>434 ± 106</td><td>235 ± 107</td><td>118 ± 81</td><td>205 ± 79</td><td>153 ± 89</td><td>180 ± 54</td><td>160 ± 86</td></tr><tr><td>IRAC 4.5 µm</td><td>598 ± 130</td><td>275 ± 102</td><td>310 ± 80</td><td>92 ± 100</td><td>363 ± 78</td><td>376 ± 67</td><td>186 ± 88</td></tr><tr><td>IRAC 5.8 µm</td><td>-643 ± 3000</td><td>488 ± 3604</td><td>-1188 ± 2445</td><td>1378 ± 2596</td><td>368 ± 2423</td><td>-699 ± 3257</td><td>1815 ± 2016</td></tr><tr><td>IRAC 8.0 µm</td><td>-3325 ± 3803</td><td>3218 ± 5706</td><td>2049 ± 4168</td><td>-4049 ± 2843</td><td>-207 ± 3771</td><td>2933 ± 3972</td><td>1665 ± 3490</td></tr></tbody></table>
<table><caption>**Table 9**
Flux Density Measurements for UVISTA-Y13 to UVISTA-Y16</caption><thead><tr><th>Filter<br>Name</th><th>UVISTA-Y13<br>(nJy)</th><th>UVISTA-Y14<br>(nJy)</th><th>UVISTA-Y15<br>(nJy)</th><th>UVISTA-Y16<br>(nJy)</th></tr></thead><tbody><tr><td>CFHTLS <i>u</i>*</td><td>-6 ± 13</td><td>-23 ± 18</td><td>...</td><td>-28 ± 18</td></tr><tr><td>SSC <i>B</i></td><td>-2 ± 6</td><td>-4 ± 7</td><td>8 ± 9</td><td>-6 ± 11</td></tr><tr><td>HSC <i>g</i></td><td>-10 ± 14</td><td>4 ± 17</td><td>-7 ± 19</td><td>13 ± 16</td></tr><tr><td>CFHTLS <i>g</i></td><td>-2 ± 10</td><td>-1 ± 15</td><td>...</td><td>-4 ± 15</td></tr><tr><td>SSC <i>V</i></td><td>-11 ± 17</td><td>11 ± 20</td><td>-5 ± 21</td><td>18 ± 24</td></tr><tr><td>HSC <i>r</i></td><td>5 ± 11</td><td>12 ± 13</td><td>-18 ± 17</td><td>4 ± 16</td></tr><tr><td>CFHTLS <i>r</i></td><td>4 ± 15</td><td>2 ± 19</td><td>...</td><td>-4 ± 22</td></tr><tr><td>SSC <i>r</i><sup>+</sup></td><td>3 ± 15</td><td>4 ± 18</td><td>-11 ± 19</td><td>-18 ± 23</td></tr><tr><td>SSC <i>i</i><sup>+</sup></td><td>3 ± 24</td><td>-1 ± 28</td><td>22 ± 29</td><td>10 ± 35</td></tr><tr><td>CFHTLS <i>y</i></td><td>14 ± 18</td><td>23 ± 26</td><td>...</td><td>-16 ± 27</td></tr><tr><td>CFHTLS <i>i</i></td><td>-1 ± 17</td><td>8 ± 24</td><td>...</td><td>1 ± 31</td></tr><tr><td>HSC <i>i</i></td><td>0 ± 18</td><td>20 ± 22</td><td>26 ± 26</td><td>-7 ± 24</td></tr><tr><td>CFHTLS <i>z</i></td><td>-17 ± 43</td><td>54 ± 58</td><td>...</td><td>-11 ± 64</td></tr><tr><td>HSC <i>z</i></td><td>14 ± 26</td><td>17 ± 27</td><td>47 ± 38</td><td>43 ± 35</td></tr><tr><td>SSC <i>z</i><sup>+</sup></td><td>-21 ± 60</td><td>29 ± 68</td><td>-41 ± 81</td><td>14 ± 75</td></tr><tr><td>HSC <i>y</i></td><td>-4 ± 69</td><td>21 ± 81</td><td>24 ± 90</td><td>-25 ± 86</td></tr><tr><td>UVISTA <i>Y</i></td><td>38 ± 45</td><td>-37 ± 50</td><td>69 ± 56</td><td>91 ± 56</td></tr><tr><td>UVISTA <i>J</i></td><td>123 ± 47</td><td>144 ± 62</td><td>189 ± 56</td><td>281 ± 64</td></tr><tr><td>UVISTA <i>H</i></td><td>166 ± 65</td><td>206 ± 75</td><td>215 ± 73</td><td>268 ± 81</td></tr><tr><td>UVISTA <i>K</i><sub>S</sub></td><td>245 ± 66</td><td>34 ± 77</td><td>106 ± 77</td><td>259 ± 93</td></tr><tr><td>IRAC 3.6 µm</td><td>177 ± 96</td><td>191 ± 69</td><td>134 ± 70</td><td>221 ± 79</td></tr><tr><td>IRAC 4.5 µm</td><td>362 ± 94</td><td>210 ± 93</td><td>-7 ± 83</td><td>389 ± 82</td></tr><tr><td>IRAC 5.8 µm</td><td>-287 ± 3140</td><td>-1202 ± 3568</td><td>1349 ± 2446</td><td>-2535 ± 2598</td></tr><tr><td>IRAC 8.0 µm</td><td>7391 ± 4246</td><td>1553 ± 4921</td><td>-1814 ± 3403</td><td>2196 ± 3371</td></tr></tbody></table>

<caption>**Table 10**</caption>

$V_{\text{max}}$ Determination of the UV LF When UVISTA-Y3 Is Treated as a Single Source

<table><thead><tr><th>$M_{\text{UV}}$<br>(mag)</th><th>$\phi$<br>(×10<sup>-3</sup> mag<sup>-1</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>-22.55</td><td>0.0008<sup>+0.0007</sup><sub>-0.0004</sub></td></tr><tr><td>-22.05</td><td>0.0017<sup>+0.0012</sup><sub>-0.0007</sub></td></tr><tr><td>-21.55</td><td>0.0036<sup>+0.0018</sup><sub>-0.0012</sub></td></tr></tbody></table>

## Appendix C

### LF Estimate When UVISTA-Y3 Is Considered as One Single Source

In this section we present, for completeness, $V_{\text{max}}$ measurements of the $z \sim 8$ UV LF when UVISTA-Y3 is considered as a single source. Table 10 lists the number densities, while in Figure 22 we present these measurements and compare them to the estimates obtained in Section 6.5 assuming a multiple component nature of UVISTA-Y3. While the impact of removing three sources from the lowest luminosity bin while adding a source to a higher luminosity bin should be obvious, our treatment of this source does not change either of the impacted LF points by more than $1\sigma$.

## Appendix D

### LF Estimate after Correcting for Lensing Magnification

Here we present the $z \sim 8$ LF determination after correcting for lensing magnification the luminosities of four objects

<figure><img src="image_25.png" /><figcaption>**Figure 22.** Comparison between the LF estimate obtained after deblending UVISTA-Y3 (filled blue circles) and that when considering UVISTA-Y3 as a single source (open symbols). Previous LF determinations at $z \sim 8$ from Bouwens et al. (2015) are presented for comparison. The systematic differences are within the $1\sigma$ uncertainties.</figcaption></figure>

(Y6, Y8, Y9, and Y13—see Section 5.5). None of the sources affected by lensing satisfies the criteria for inclusion in the $z \sim 9$ LF.

The change of solid angle introduced by the lenses affects the volume estimates and therefore the volume densities. A proper estimate of the volumes that takes into account the lensing effects would require running a simulation with the actual lensing effects from all the sources in the region considered for our search. While this goes beyond the scope of this paper, we can attempt to estimate the average effect based
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Comparison between the LF estimate obtained after correcting sources for lens magnification (open orange circles) to those without any corrections (filled blue circles). Error bars for the LF without lensing correction were omitted for clarity. Previous LF determinations at $z \sim 8$ from Bouwens et al. (2015) are presented for comparison. The systematic differences are well within the $1\sigma$ uncertainties.</figcaption></figure>

<table><caption>**Table 11**
$V_{\text{max}}$ Determination of the UV LF When Accounting for Lensing Magnification</caption><thead><tr><th>$M_{\text{UV}}$<br>(mag)</th><th>$\phi$<br>(×10<sup>-3</sup> mag<sup>-1</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>-22.55</td><td>0.0008<sup>+0.0007</sup><sub>-0.0004</sub></td></tr><tr><td>-22.05</td><td>0.0010<sup>+0.0010</sup><sub>-0.0006</sub></td></tr><tr><td>-21.55</td><td>0.0037<sup>+0.0019</sup><sub>-0.0013</sub></td></tr></tbody></table>

on the statistics in our sample. Considering that 4/18 objects in our sample have cumulative lensing magnification factors of ~1.85, 1.96, 1.39, and 1.60 respectively, conserving the surface brightness would imply shrinking by the same amounts the area corresponding to each source (i.e., ~1/18 the area of the UltraVISTA field). Given that the remaining 14/18 sources are not significantly magnified, the total area available for the search would then be 91% of the area when no lensing is considered. The volume densities would then be ~10% higher (just ~0.04 dex) than those computed without introducing the lensing magnification. Because of the very marginal effect of this correction compared to the Poissonian and cosmic variance uncertainties, the LF measurements presented in this section were computed with the same volumes adopted for the LFs presented in Section 6.5, which do not include any correction for lensing effects. In Table 11 we list the corresponding $V_{\text{max}}$ determination when the luminosities of those sources are corrected for lensing magnification.

In Figure 23 we compare the $V_{\text{max}}$ estimates of the delensed sample to those from the original sample and the corresponding Schechter parameterization. A fit to the Schechter form gives $M^* = -20.81_{-0.30}^{+0.27}$ mag, $\alpha = -2.09_{-0.19}^{+0.20}$, and $\log \Phi^* = -3.84_{-0.32}^{+0.26}$. The 68% and 95% confidence intervals of the Schechter parameters are presented in Figure 24. The LF determination obtained accounting for lensing magnification is consistent (at $1\sigma$) with that obtained without such a correction. This suggests that, at least for the small samples that are the subject of the current study, a full accounting for the lensing magnification does not significantly impact on our conclusions.

<figure><img src="image_27.png" /><figcaption>**Figure 24.** 68% and 95% likelihood contours on the Schechter parameters at $z \sim 8$ derived in the present work after correcting the luminosities of the sample for lensing magnification (orange curves) compared to those from the analysis without lensing magnification correction (blue curves).</figcaption></figure>
ORCID iDs

Mauro Stefanon https://orcid.org/0000-0001-7768-5309

Rychard J. Bouwens https://orcid.org/0000-0002-4989-2471

Pascal Oesch https://orcid.org/0000-0001-5851-6649

Matthew L. N. Ashby https://orcid.org/0000-0002-3993-0745

Karina I. Caputi https://orcid.org/0000-0001-8183-1460

Marijn Franx https://orcid.org/0000-0002-8871-3026

Johan P. U. Fynbo https://orcid.org/0000-0002-8149-8298

Danilo Marchesini https://orcid.org/0000-0001-9002-3502

Henry J. McCracken https://orcid.org/0000-0002-9489-7765

Bo Milvang-Jensen https://orcid.org/0000-0002-2281-2785

Adam Muzzin https://orcid.org/0000-0002-9330-9108

Pieter van Dokkum https://orcid.org/0000-0002-8282-9888

References

Aihara, H., Arimoto, N., Armstrong, R., et al. 2018b, PASJ, 70, 4

Aihara, H., Armstrong, R., Bickerton, S., et al. 2018a, PASJ, 70, 8

Anders, P., & Fritze-v Alvensleben, U. 2003, A&A, 401, 1063

Ashby, M. L. N., Caputi, K. I., Cowley, W., et al. 2018, ApJS, 237, 39

Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2013, ApJ, 769, 80

Ashby, M. L. N., Willner, S. P., Fazio, G. G., et al. 2015, ApJS, 218, 33

Atek, H., Siana, B., Scarlata, C., et al. 2011, ApJ, 743, 121

Avni, Y., & Bahcall, J. N. 1980, ApJ, 235, 694

Barone-Nugent, R. L., Wyithe, J. S. B., Trenti, M., et al. 2015, MNRAS, 450, 1224

Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729

Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212

Bernard, S. R., Carrasco, D., Trenti, M., et al. 2016, ApJ, 827, 76

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Bolton, J. S., & Haehnelt, M. G. 2007, MNRAS, 374, 493

Bouwens, R. J., Bradley, L., Zitrin, A., et al. 2014, ApJ, 795, 126

Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2008, ApJ, 686, 230

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2010, ApJL, 708, L69

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2011, ApJ, 737, 90

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bouwens, R. J., Oesch, P. A., Labbé, I., et al. 2016, ApJ, 830, 67

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., & McLeod, D. J. 2017, MNRAS, 466, 3612

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2012, MNRAS, 426, 2772

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2014, MNRAS, 440, 2810

Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2015, MNRAS, 452, 1817

Bradley, L. D., Trenti, M., Oesch, P. A., et al. 2012, ApJ, 760, 108

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bridge, J. S., Holwerda, B. W., Stefanon, M., et al. 2019, arXiv:1907.05512

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Burgasser, A. J. 2014, ASINc, 11, 7

Calvi, V., Trenti, M., Stiavelli, M., et al. 2016, ApJ, 817, 120

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Caputi, K. I., Deshmukh, S., Ashby, M. L. N., et al. 2017, ApJ, 849, 45

Castellano, M., Fontana, A., Grazian, A., et al. 2012, A&A, 540, A39

Chabrier, G. 2003, PASP, 115, 763

de Barros, S., Oesch, P. A., Labbé, I., et al. 2018, MNRAS, 489, 2355

Eldridge, J. J., Stanway, E. R., Xiao, L., et al. 2017, PASA, 34, e058

Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJL, 763, L7

Erb, D. K., Pettini, M., Shapley, A. E., et al. 2010, ApJ, 719, 1168

Erben, T., Hildebrandt, H., Lerchster, M., et al. 2009, A&A, 493, 1197

Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122

Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10

Fialkov, A., & Loeb, A. 2015, ApJ, 806, 256

Finkelstein, S. L., Papovich, C., Dickinson, M., et al. 2013, Natur, 502, 524

Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164

Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71

Finlator, K., Oppenheimer, B. D., & Davé, R. 2011, MNRAS, 410, 1703

Galametz, A., Grazian, A., Fontana, A., et al. 2013, ApJS, 206, 10

Guo, Y., Ferguson, H. C., Giavalisco, M., et al. 2013, ApJS, 207, 24

Harikane, Y., Ouchi, M., Ono, Y., et al. 2018, PASJ, 70, S11

Hildebrandt, H., Pielorz, J., Erben, T., et al. 2009, A&A, 498, 725

Hoag, A., Bradač, M., Brammer, G., et al. 2018, ApJ, 854, 39

Hoag, A., Bradač, M., Trenti, M., et al. 2017, NatAs, 1, 0091

Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6

Huang, K.-H., Lemaux, B. C., Schmidt, K. B., et al. 2016, ApJL, 823, L14

Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6

Inoue, A. K. 2011, MNRAS, 415, 2920

Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2018, ApJ, 854, 73

Izotov, Y. I., Chaffee, F. H., Foltz, C. B., et al. 1999, ApJ, 527, 757

Jaacks, J., Nagamine, K., & Choi, J. H. 2012, MNRAS, 427, 403

Kochanek, C. S., Schneider, P., & Wambsganss, J. 2004, arXiv:astro-ph/0407232

Kormann, R., Schneider, P., & Bartelmann, M. 1994, A&A, 284, 285

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Kron, R. G. 1980, ApJS, 43, 305

Labbé, I., Bouwens, R., Illingworth, G. D., & Franx, M. 2006, ApJL, 649, L67

Labbé, I., González, V., Bouwens, R. J., et al. 2010a, ApJL, 716, L103

Labbé, I., González, V., Bouwens, R. J., et al. 2010b, ApJL, 708, L26

Labbé, I., Oesch, P. A., Bouwens, R. J., et al. 2013, ApJL, 777, L19

Labbé, I., Oesch, P. A., Illingworth, G. D., et al. 2015, ApJS, 221, 23

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Laporte, N., Streblyanska, A., Clement, B., et al. 2014, A&A, 562, L8

Laporte, N., Streblyanska, A., Kim, S., et al. 2015, A&A, 575, A92

Larson, R. L., Finkelstein, S. L., Pirzkal, N., et al. 2018, ApJ, 858, 94

Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219

Livermore, R. C., Trenti, M., Bradley, L. D., et al. 2018, ApJL, 861, L17

Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97

Marsan, Z. C., Marchesini, D., Muzzin, A., et al. 2019, ApJ, 871, 201

Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101

Mason, C. A., Trenti, M., & Treu, T. 2015, ApJ, 813, 21

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McGuire, J. T. W., Tanvir, N. R., Levan, A. J., et al. 2016, ApJ, 825, 135

McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812

McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2015, MNRAS, 450, 3032

McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Morishita, T., Trenti, M., Stiavelli, M., et al. 2018, ApJ, 867, 150

Moster, B. P., Somerville, R. S., Newman, J. A., & Rix, H.-W. 2011, ApJ, 731, 113

Mowla, L., van Dokkum, P., Brammer, G., et al. 2019, ApJ, 880, 57

Murray, S. G., Power, C., & Robotham, A. S. G. 2013, A&C, 3, 75

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJS, 206, 8

Nayyeri, H., Hemmati, S., Mobasher, B., et al. 2017, ApJS, 228, 7

Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010, ApJ, 709, L21

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2012, ApJ, 759, 135

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2014, ApJ, 786, 108

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015a, ApJ, 808, 104

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105

Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129

Oesch, P. A., Carollo, C. M., Stiavelli, M., et al. 2009, ApJ, 690, 1350

Oesch, P. A., Labbé, I., Bouwens, R. J., et al. 2013, ApJ, 772, 136

Oesch, P. A., van Dokkum, P. G., Illingworth, G. D., et al. 2015b, ApJL, 804, L30

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Ono, Y., Ouchi, M., Curtis-Lake, E., et al. 2013, ApJ, 777, 155

Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10

Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, ApJ, 744, 83

Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, ApJ, 743, 132

Pentericci, L., Vanzella, E., Castellano, M., et al. 2018, A&A, 619, A147

Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, ApJ, 793, 113

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Rasappu, N., Smit, R., Labbé, I., et al. 2016, MNRAS, 461, 3886

Roberts-Borsani, G. W., Bouwens, R. J., Oesch, P. A., et al. 2016, ApJ, 823, 143

Robertson, B. E., Ellis, R. S., Dunlop, J. S., McLure, R. J., & Stark, D. P. 2010, Natur, 468, 49

Rogers, A. B., McLure, R. J., Dunlop, J. S., et al. 2014, MNRAS, 440, 3714

Salmon, B., Coe, D., Bradley, L., et al. 2017, arXiv:1710.08930

Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Schechter, P. 1976, ApJ, 203, 297
Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, ApJ, 795, 20

Schenker, M. A., Robertson, B. E., Ellis, R. S., et al. 2013, ApJ, 768, 196

Schenker, M. A., Stark, D. P., Ellis, R. S., et al. 2012, ApJ, 744, 179

Schmidt, K. B., Treu, T., Bradač, M., et al. 2016, ApJ, 818, 38

Schmidt, K. B., Treu, T., Trenti, M., et al. 2014, ApJ, 786, 57

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1

Shibuya, T., Kashikawa, N., Ota, K., et al. 2012, ApJ, 752, 114

Shibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15

Shim, H., Chary, R.-R., Dickinson, M., et al. 2011, ApJ, 738, 69

Shull, J. M., Harness, A., Trenti, M., & Smith, B. D. 2012, ApJ, 747, 100

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Smit, R., Bouwens, R. J., Carniani, S., et al. 2018, Natur, 553, 178

Smit, R., Bouwens, R. J., Franx, M., et al. 2015, ApJ, 801, 122

Smit, R., Bouwens, R. J., Labbé, I., et al. 2014, ApJ, 784, 58

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

Stark, D. P. 2016, ARA&A, 54, 761

Stark, D. P., Ellis, R. S., Charlot, S., et al. 2017, MNRAS, 464, 469

Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, ApJL, 728, L2

Stark, D. P., Richard, J., Charlot, S., et al. 2015b, MNRAS, 450, 1846

Stark, D. P., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129

Stark, D. P., Walth, G., Charlot, S., et al. 2015a, MNRAS, 454, 1393

Stefanon, M., Labbé, I., Bouwens, R. J., et al. 2017b, ApJ, 851, 43

Stefanon, M., Yan, H., Mobasher, B., et al. 2017a, ApJS, 229, 32

Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68

Szomoru, D., Franx, M., van Dokkum, P. G., et al. 2013, ApJ, 763, 73

Taniguchi, Y., Scoville, N., Murayama, T., et al. 2007, ApJS, 172, 9

Tanvir, N. R., Fynbo, J. P. U., de Ugarte Postigo, A., et al. 2019, MNRAS, 483, 5380

Tanvir, N. R., Levan, A. J., Fruchter, A. S., et al. 2012, ApJ, 754, 46

Tilvi, V., Papovich, C., Finkelstein, S. L., et al. 2014, ApJ, 794, 5

Trac, H., Cen, R., & Mansfield, P. 2015, ApJ, 813, 54

Trenti, M., Bradley, L. D., Stiavelli, M., et al. 2011, ApJL, 727, L39

Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013, ApJL, 775, L29

van der Wel, A., Chang, Y.-Y., Bell, E. F., et al. 2014, ApJL, 792, L6

van der Wel, A., Straughn, A. N., Rix, H.-W., et al. 2011, ApJ, 742, 111

Vanzella, E., Pentericci, L., Fontana, A., et al. 2011, ApJL, 730, L35

Waters, D., Wilkins, S. M., Di Matteo, T., et al. 2016, MNRAS, 461, L51

Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. 2011, Natur, 469, 181

Yan, H., Yan, L., Zamojski, M. A., et al. 2011, ApJL, 728, L22

Zitrin, A., Labbé, I., Belli, S., et al. 2015, ApJL, 810, L12