# HST Imaging of the Brightest z ~ 8-9 Galaxies from UltraVISTA: The Extreme Bright End of the UV Luminosity Function

Mauro Stefanon, Ivo Labbé, Rychard J. Bouwens, Gabriel B. Brammer, Pascal Oesch, Marijn Franx, Johan P. U. Fynbo, Bo Milvang-Jensen, Adam Muzzin, Garth D. Illingworth, et al.

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

To cite this version:

Mauro Stefanon, Ivo Labbé, Rychard J. Bouwens, Gabriel B. Brammer, Pascal Oesch, et al. HST Imaging of the Brightest z ~ 8-9 Galaxies from UltraVISTA: The Extreme Bright End of the UV Luminosity Function. The Astrophysical Journal, 2017, 851, 10.3847/1538-4357/aa9a40. insu-03666142

HAL Id: insu-03666142

https://insu.hal.science/insu-03666142v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# HST Imaging of the Brightest $z \sim 8-9$ Galaxies from UltraVISTA: The Extreme Bright End of the UV Luminosity Function

Mauro Stefanon¹, Ivo Labbé¹, Rychard J. Bouwens¹, Gabriel B. Brammer², Pascal Oesch³, Marijn Franx¹,
Johan P. U. Fynbo⁴, Bo Milvang-Jensen⁴, Adam Muzzin⁵, Garth D. Illingworth⁶, Olivier Le Fèvre⁷, Karina I. Caputi⁸,
Benne W. Holwerda⁹, Henry J. McCracken¹⁰, Renske Smit¹¹,¹², and Dan Magee⁶

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

²Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

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

⁴Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100 Copenhagen Ø, Denmark

⁵York University, 4700 Keele Street, Toronto, ON, M3J 1P3, Canada

⁶UCO/Lick Observatory, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA

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

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

⁹Department of Physics and Astronomy, University of Louisville, Louisville KY 40292, USA

¹⁰Institut d'Astrophysique de Paris, 98bis Boulevard Arago, F-75014 Paris, France

¹¹Cavendish Laboratory, University of Cambridge, 19 JJ Thomson Avenue, Cambridge CB3 0HE, UK

¹²Kavli Institute of Cosmology c/o Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK

Received 2017 June 14; revised 2017 November 8; accepted 2017 November 9; published 2017 December 11

## Abstract

We report on the discovery of three especially bright candidate $z_{phot} \gtrsim 8$ galaxies. Five sources were targeted for follow-up with the *Hubble Space Telescope (HST)/Wide Field Camera 3 (WFC3)*, selected from a larger sample of 16 bright ($24.8 \lesssim H \lesssim 25.5$ mag) candidate $z \gtrsim 8$ Lyman break galaxies (LBGs) identified over 1.6 degrees² of the COSMOS/UltraVISTA field. These were selected as *Y* and *J* dropouts by leveraging the deep (*Y*-to-*Ks* ~ 25.3–24.8 mag, 5σ) NIR data from the UltraVISTA DR3 release, deep ground-based optical imaging from the CFHTLS and Suprime-Cam programs, and *Spitzer/IRAC* mosaics combining observations from the SMUVS and SPLASH programs. Through the refined spectral energy distributions, which now also include new HyperSuprimeCam *g*, *r*, *i*, *z*, and *Y*-band data, we confirm that 3/5 galaxies have robust $z_{phot} \sim 8.0–8.7$, consistent with the initial selection. The remaining 2/5 galaxies have a nominal $z_{phot} \sim 2$. However, with *HST* data alone, these objects have increased probability of being at $z \sim 9$. We measure mean UV continuum slopes $\beta = -1.74 \pm 0.35$ for the three $z \sim 8-9$ galaxies, marginally bluer than similarly luminous $z \sim 4-6$ in CANDELS but consistent with previous measurements of similarly luminous galaxies at $z \sim 7$. The circularized effective radius for our brightest source is $0.9 \pm 0.3$ kpc, similar to previous measurements for a bright $z \sim 11$ galaxy and bright $z \sim 7$ galaxies. Finally, enlarging our sample to include the six brightest $z \sim 8$ LBGs identified over UltraVISTA (i.e., including three other sources from Labbé et al.) we estimate for the first time the volume density of galaxies at the extreme bright end ($M_{UV} \sim -22$ mag) of the $z \sim 8$ UV luminosity function. Despite this exceptional result, the still large statistical uncertainties do not allow us to discriminate between a Schechter and a double-power-law form.

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

## 1. Introduction

The study of the galaxy populations at the epoch of re-ionization has substantially improved in the last decade thanks to the exceptional sensitivity of the *Hubble Space Telescope (HST)/Wide Field Camera 3 (WFC3)*. Programs such as the Hubble Deep Field (Williams et al. 1996), the Hubble Ultra-Deep and eXtreme-Deep Field (Beckwith et al. 2006; Illingworth et al. 2013), the Brightest of the Reionizing Galaxies (BoRG, Trenti et al. 2011), the Hubble Frontier Fields (Lotz et al. 2017), and the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CANDELS)/3D-HST (Grogin et al. 2011; Koekemoer et al. 2011; van Dokkum et al. 2011; Brammer et al. 2012; Momcheva et al. 2016) have enabled the identification of ~1500 candidate galaxies at $z > 6$, ~200 of which are at $z \sim 8–10$ (e.g., Oesch et al. 2013; Bouwens et al. 2015, 2016; Finkelstein et al. 2015). These samples are characterized by $M_{UV} \gtrsim -22$ (i.e., apparent magnitudes fainter than $\gtrsim 25.5$ at $z \sim 8$), ~1 mag more

luminous than the current determinations of the characteristic magnitude of the rest-frame UV luminosity functions (LFs).

Given the steep faint-end slope of the UV LF at $z \gtrsim 6$ (Schechter 1976) ($\alpha \sim -2$; see e.g., Bouwens et al. 2011, 2015; McLure et al. 2013; Schenker et al. 2013; Duncan et al. 2014; Finkelstein et al. 2015), galaxies fainter than the characteristic luminosity dominate the estimates of the star formation rate (SFR) density (e.g., Oesch et al. 2014; McLeod et al. 2015); furthermore, under the hypothesis that the faint-end slope does not decrease at luminosities 3–4 mag fainter than current observational limits at $z \sim 6–8$, their much higher (factors of $10^2–10^4$) volume density, compared to the bright end, has been proven sufficient for low-luminosity galaxies to complete the re-ionization by $z \sim 6$ (e.g., Stark 2016 and references therein). Nonetheless, given the correlation between the faint-end slope and the characteristic magnitude of the Schechter parameterization, commonly adopted to describe the shape of the LF at high redshift, the determination of the faint-
end slope also benefits from improvements at the bright side. Furthermore, the identification of luminous galaxies at early epochs constitutes an important constraint to all models of galaxy formation and evolution. The recent spectroscopic confirmation of GN-z11, a luminous ($M_{UV} = -22.1$ mag) galaxy at the record redshift of $z_{grism} = 11.09^{+0.08}_{-0.12}$, has shown that its associated number density is higher than both the extrapolation to $z \sim 10$ of the Schechter parameterization of the UV LFs and the model predictions (Oesch et al. 2016), challenging our current understanding of galaxy formation and evolution.

The steep exponential decline at the bright end of the current UV LF determinations suggests that probing the LF at even brighter luminosities requires exploring areas of the order of a square degree in NIR bands to depths of $\sim$25 mag. Some progress in this direction has come from the BoRG and HIPPIES programs (Trenti et al. 2011; Yan et al. 2011; Bernard et al. 2016; Calvi et al. 2016), which have uncovered galaxies up to $z \sim 8-10$ with $M_{UV} \sim -22.5$ mag (e.g., Calvi et al. 2016). A complementary approach comes from ground-based surveys, which allow us to extend the surveyed area to $\sim$1–100 deg$^2$, necessary to minimize the effects of cosmic variance in the systematic search for the brightest objects.

Recently, Bowler et al. (2014, 2015, 2017) identified a sample of luminous galaxies ($-23 \lesssim M_{UV} \lesssim -22$ mag) at $z \sim 6-7$ in the COSMOS/UltraVISTA field. Interestingly, the associated number densities are in excess of the Schechter (1976) form, suggesting that a double power law could provide a better description at the bright end than the commonly assumed Schechter form. Partial confirmation of this comes from Ono et al. (2017) who measured the bright end of the UV LF at $4 \lesssim z \lesssim 7$ using data from the HyperSuprimeCam Survey (Aihara et al. 2017a, 2017b). Analysis of this three-layered data set resulted in a sample of $\sim 600 z \sim 6-7$ Lyman break galaxies (LBGs) ($\sim 70$ galaxies at $z \sim 7$) with $M_{UV} \lesssim -25$ mag over $\sim 100$ deg$^2$. After carefully removing active galactic nucleus (AGN) contaminants, their UV LF measurements show an excess at the bright end of the UV LF compared to the Schechter parameterization from previous studies, although a double power law still over-predicts the brightest end.

In order to probe the bright end of the UV LF at even higher redshift we leveraged the deep and ultradeep data of the third data release (DR3) of the UltraVISTA program (McCracken et al. 2012), complemented by deep optical data from the Canada-France-Hawaii Legacy Survey (CFHTLS; Erben et al. 2009; Hildebrandt et al. 2009) and Subaru/Suprime-Cam (SSC; Taniguchi et al. 2007), and with deep IRAC mosaics that we generated following Labbé et al. (2015) which include observations from the SMUVS (PI: K. Caputi) and SPLASH (PI: P. Capak) programs. Using LBG criteria we selected a sample of 16 bright ($H \sim 24-25$ AB) galaxies at $z \sim 8$ (I. Labbé et al. 2017, in preparation).

The primary question is whether the bright sources identified from the ground-based selections exist or whether they are a population of lower-$z$ interlopers. Indeed, spectroscopic confirmation has recently been obtained for three UV-luminous ($M_{UV} \sim -22$ mag) galaxies at $z \sim 7.5-8.7$ with $H \sim 25.1$ AB selected from CANDELS fields (Roberts-Borsani et al. 2016), one at $z \sim 8.7$ (Zitrin et al. 2015) and one at $z \sim 11$ (GN-z11; Oesch et al. 2016). Furthermore, the lower spatial resolution of ground-based observations, compared to HST data, can blend

the signal from mergers or from physically unrelated objects and hence make them appear as single sources (e.g., Bowler et al. 2017), resulting in an over-estimate of the bright end of the LF and an underestimate of the LF at lower luminosities. Photometric variability can be indicative of the presence of an AGN component or of a stellar or brown dwarf contaminant, which would introduce systematics or even contaminate the sample. Fluctuations in the signal induced by the random noise from the background can potentially conspire to suppressing low signal-to-noise ratio (S/N) signal at optical wavelengths and thus mimicking a high-redshift solution. Moreover, Bowler et al. (2017) have shown that the electronics of the detectors can introduce image ghosts that can mimic high-redshift objects. While each of these effects are likely rare, we are looking for a small number of real high-redshift candidates in a large imaging data set, and follow-up imaging is required to validate these candidates, effectively eliminating many of these concerns.

We therefore selected five of the brightest candidate $z \sim 8-9$ LBGs from I. Labbé et al. (2017, in preparation) for targeted follow-up with HST/WFC3 $z_{098}$, $J_{125}$, and $H_{160}$ bands (PI: R. Bouwens, PID: 14895) to attempt to confirm these sources and to better constrain their physical properties. These candidates stood out for their unprecedented brightness ($24.5 \lesssim H \lesssim 25.2$) and for their tantalizing plausible $z_{phot} \sim 8.5-9.0$ solution, being detected in the UltraVISTA ultradeep stripes and non-detection in the deepest optical ground-based data.

This paper is devoted to presenting the results of the analysis of the HST data for the five candidate $z \sim 8$ galaxies. In Section 2 we briefly describe the data sets and the selection criteria adopted for the assembly of the sample; in Section 3 we describe the HST data and how the photometry was performed; the results are presented in Section 4; in Section 5 we discuss the results and in Section 6 we conclude.

Throughout this work, we use the following short form to indicate HST/WFC3 filters: WFC3/F098M → $z_{098}$; WFC3/F105W → $Y_{105}$; WFC3/F125W → $J_{125}$; WFC3/F140W → $JH_{140}$; and WFC3/F160W → $H_{160}$. We adopt a cosmology with $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_{\Lambda} = 0.7$, and $\Omega_m = 0.3$. All magnitudes are in the AB system.

## 2. Sample Selection

In this section we briefly summarize the data set and the procedure followed to select the sample of candidate luminous $z \sim 8$ sources in the COSMOS/UltraVISTA field. Full details are given in an accompanying paper (I. Labbé et al. 2017, in preparation). We give a concise summary below.

Our sample is based on the deep NIR imaging available over the COSMOS field (Scoville et al. 2007) from the third release (DR3)<fnref n="13" /> of the UltraVISTA program (McCracken et al. 2012). This data release provides mosaics in the *Y*, *J*, *H*, and *K*<sub>S</sub> broad bands together with a narrowband centered at 1.18 µm (NB118). The mosaics in the broadband filters are characterized by four ultradeep stripes reaching *Y*-to-*K*<sub>S</sub> $\sim 24.8-25.3$ mag ($5\sigma$, 1"/2 aperture diameter corrected to total), alternating with four deep stripes (*Y*-to-*K*<sub>S</sub> $\sim 23.7-24.4$ mag, $5\sigma$, 1"/2 aperture to total), for a total area of $\sim 1.6$ degree². The UltraVISTA data were complemented by deep optical imaging from CFHT/Megacam in *g*, *r*, *i* and *z* (Erben et al. 2009; Hildebrandt et al. 2009) from

<footnote n="13">¹³ https://www.eso.org/sci/observing/phase3/data_releases/uvista_dr3.pdf</footnote>
<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.3</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>CFHTLS <i>g</i></td><td>2.2</td><td>26.8</td></tr><tr><td>SSC <i>V</i></td><td>2.2</td><td>26.4</td></tr><tr><td>HSC <i>r</i><sup>c</sup></td><td>1.6</td><td>26.8</td></tr><tr><td>CFHTLS <i>r</i></td><td>2.1</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>CFHTLS <i>y</i></td><td>2.0</td><td>26.1</td></tr><tr><td>CFHTLS <i>i</i></td><td>2.0</td><td>26.0</td></tr><tr><td>HSC <i>i</i><sup>c</sup></td><td>1.7</td><td>26.3</td></tr><tr><td>CFHTLS <i>z</i></td><td>2.1</td><td>25.2</td></tr><tr><td>HSC <i>z</i><sup>c</sup></td><td>1.6</td><td>25.9</td></tr><tr><td>SSC <i>z</i><sup>+</sup></td><td>2.3</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.6</td><td>25.4/24.5</td></tr><tr><td>UVISTA <i>J</i></td><td>2.4</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.2</td><td>24.8/23.7</td></tr><tr><td>IRAC 3.6 µm</td><td>5.3</td><td>24.9/24.5</td></tr><tr><td>IRAC 4.5 µm</td><td>5.4</td><td>24.7/24.3</td></tr><tr><td>IRAC 5.8 µm</td><td>8.4</td><td>20.8</td></tr><tr><td>IRAC 8.0 µm</td><td>10.1</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 two 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 HSC 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.

the CFHTLS and SSC in $B_j$, $V_j$, $r^+$, $i^+$, and $z$ bands (Taniguchi et al. 2007). Full-depth mosaics were constructed following Labbé et al. (2015) for the *Spitzer*/IRAC 3.6 and 4.5 µm observations from S-COSMOS (Sanders et al. 2007), the *Spitzer* Extended Deep Survey (SEDS; Ashby et al. 2013), the *Spitzer* Cosmic Assembly Near-Infrared Deep Extragalactic Survey (S-CANDELS, Ashby et al. 2015), *Spitzer* Large Area Survey with HyperSuprimeCam (HSC) (SPLASH, PI: Capak), and the complete set of observations of the *Spitzer* Matching Survey of the UltraVISTA ultradeep Stripes (SMUVS, PI: Caputi; M. L. N. Ashby et al. 2017, in preparation; Caputi et al. 2017).

Table 1 lists the 5σ depths of the adopted ground-based and *Spitzer*/IRAC mosaics. They were measured as the standard deviation of fluxes in ~4000 empty apertures of 1"/2 diameter, randomly scattered across the mosaic avoiding sources in the segmentation map. The values were finally multiplied by the average aperture corrections for each band (also reported in Table 1) to convert them into total fluxes, mimicking the procedures adopted for the flux measurements. While the exposure times across the CFHTLS, SSC, and HSC mosaics are fairly uniform, the UltraVISTA and IRAC 3.6 and 4.5 µm

mosaics are roughly characterized by a bi-modal depth. In these bands we therefore computed two different depths, restricting the random locations to regions representative of either one of the two typical depths. Our depth measurements are ≈0.5–0.8 mag brighter than previous estimates (e.g., Bowler et al. 2014; Skelton et al. 2014). One possible reason for this is the specific statistical estimators adopted for the measurement. For instance, Bowler et al. (2014) compute the background noise using the median absolute deviation (MAD). For a normal distribution, MAD is a factor ~1.5 lower than the standard deviation, thus corresponding to ~0.4 mag fainter estimates. Finally, to ensure basic consistency with the results of Bowler et al. (2014), we independently made use of the MAD estimator to measure 5σ depths and recovered values within 0.1 mag from those presented by Bowler et al. (2014).<fnref n="14" />

Our search was carried out on the whole 1.6 degree² of the UltraVISTA field. The mosaics of the optical and ground-based NIR bands were point-spread function (PSF) homogenized to the UltraVISTA *J* band, so that the flux curve of growth for a point source would be the same across all bands. Fluxes from these bands were extracted using SEExtractor (Bertin & Arnouts 1996) in dual mode. Source detection was performed on the square root of the χ² image (Szalay et al. 1999) built using the UltraVISTA *J*-, *H*-, and *K*<sub>S</sub>-band science and rms-map mosaics. Total fluxes were computed from 1"/2-diameter apertures and applying a correction based on the PSF and brightness profile of each individual object.

Flux measurement for the *Spitzer*/IRAC bands was performed with the *mophongo* software (Labbé et al. 2006, 2010a, 2010b, 2013, 2015); briefly, the procedure consists in reconstructing the light profile of the objects in the same field of the source under consideration, using as a prior the morphological information from a higher-resolution image. Successively, all the neighboring objects within a radius of 15"/0 from the source under analysis are removed using the positional and morphological information from the high-resolution image and a careful reconstruction of the convolution kernel (see also e.g., Fernández-Soto et al. 1999; Laidler et al. 2007; Merlin et al. 2015). Finally, aperture photometry is performed on the neighbor-cleaned source. For this work, we adopted an aperture of 1"/2 diameter. The model profile of the individual sources is finally used to correct the aperture fluxes for missing light outside the aperture. Specifically, this correction to *total* flux is performed irrespective of detections or non-detections/negative flux measurements.

We note here that the use of morphological information and the kernel reconstruction operated by *mophongo* (similarly to other codes based on template fitting) renders unnecessary matching the images to the broadest PSF in the sample prior to extracting the flux densities, further reducing potential contaminations from neighboring sources.

The sample of candidate galaxies at $z \gtrsim 8$ was selected applying Lyman break criteria. Specifically, the following two color criteria were applied (I. Labbé et al. 2017, in preparation):

$$Y - (J + H)/2 > 0.75 \lor (J - H) > 0.8 \quad (1)$$

<footnote n="14">¹⁴ This test was performed on 1"/8 empty apertures for full consistency with Bowler et al. (2014).</footnote>
<table><caption>**Table 2**
HST Observations We Obtained over the Bright $z \sim 8-9$ Candidates from I. Labbé et al. (2017, in preparation)</caption><thead><tr><th rowspan="2">ID</th><th rowspan="2">R.A.<br>(J2000)</th><th rowspan="2">Decl.<br>(J2000)</th><th rowspan="2">H<br>(mag)</th><th colspan="5">Exposure Times</th><th colspan="5">Photometric Depths</th></tr><tr><th>z<sub>098</sub><br>(s)</th><th>Y<sub>105</sub><sup>a</sup><br>(s)</th><th>J<sub>125</sub><br>(s)</th><th>JH<sub>140</sub><sup>a</sup><br>(s)</th><th>H<sub>160</sub><br>(s)</th><th>z<sub>098</sub><br>(mag)</th><th>Y<sub>105</sub><sup>a</sup><br>(mag)</th><th>J<sub>125</sub><br>(mag)</th><th>JH<sub>140</sub><sup>a</sup><br>(mag)</th><th>H<sub>160</sub><br>(mag)</th></tr></thead><tbody><tr><td>UVISTA-Y-1</td><td>09:57:47.90</td><td>2:20:43.7</td><td>24.8</td><td>1512</td><td>1022</td><td>462</td><td>1197</td><td>412</td><td>25.8</td><td>25.6</td><td>25.0</td><td>25.5</td><td>24.6</td></tr><tr><td>UVISTA-Y-5</td><td>10:00:31.89</td><td>1:57:50.2</td><td>24.9</td><td>1512</td><td>...</td><td>462</td><td>...</td><td>412</td><td>25.8</td><td>...</td><td>25.0</td><td>...</td><td>24.6</td></tr><tr><td>UVISTA-Y-6</td><td>10:00:12.51</td><td>2:03:00.5</td><td>25.3</td><td>1512</td><td>...</td><td>462</td><td>...</td><td>412</td><td>25.8</td><td>...</td><td>25.0</td><td>...</td><td>24.6</td></tr><tr><td>UVISTA-J-1</td><td>10:02:25.48</td><td>2:29:13.6</td><td>24.6</td><td>1512</td><td>...</td><td>462</td><td>...</td><td>412</td><td>25.8</td><td>...</td><td>25.0</td><td>...</td><td>24.6</td></tr><tr><td>UVISTA-J-2</td><td>09:59:07.19</td><td>1:56:54.0</td><td>24.4</td><td>1512</td><td>...</td><td>462</td><td>...</td><td>412</td><td>25.8</td><td>...</td><td>25.0</td><td>...</td><td>24.6</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, consistent with the flux measurements in the WFC3 bands used in this work.

<sup>a</sup> Fortuitously, observations in the $Y_{105}$ and $JH_{140}$ bands are available over one $z \sim 8-9$ candidate in our program as part of the separate HST program SUSHI (PI: Nao Suzuki, PID: 14808).

for the Lyman break, and

$$ (H - K_S < 0.7) \land (K_S - [3.6] < 1.75 \lor H - [3.6] < 1.75) \qquad (2) $$

to reject sources with a red continuum red-ward of the J band, likely the result of a lower redshift dusty interloper. The symbols $\land$ and $\lor$ correspond to the logical AND and OR, respectively. Furthermore, sources showing $\ge 2\sigma$ detection in any of the ground-based data blue-ward of the Lyman break were removed from the sample. We note here that Equation (1) includes two different Lyman break criteria: the first one selects galaxies whose Lyman break enters the Y band, i.e., whose redshift is $\ge 7.5$, and the second one selects galaxies whose Lyman break enters the J band, i.e., when the redshift is $\ge 9.5$.

The sample was finally cleaned from potential brown dwarf contaminants. To this aim, we opted for not adopting SExtractor `class_star` parameter because the classification becomes uncertain at low S/N (e.g., Bertin & Arnouts 1996). To overcome this, other star/galaxy separation criteria based on SExtractor have been developed (see e.g., Holwerda et al. 2014). One of the most reliable is the effective radius versus magnitude (Ryan et al. 2011). However, in order to separate stars from galaxies, this method still requires to be applied to sources about 1.5 mag above the photometric limit. Therefore, candidate brown dwarves were identified by fitting the observed SEDs with stellar templates from the SpecX prism library (Burgasser 2014) and from Burrows et al. (2006) (which provide coverage up to $\sim 15 \mu m$ for L and T dwarves) and excluding sources with $\chi^2$ from the stellar template set lower than from the galaxy templates. The above selection criteria resulted in 16 candidate $z \ge 8$ galaxies brighter than $H = 25.8$ mag.

Out of the 16 candidate galaxies at $z \ge 8$, we selected five (labeled UVISTA-Y-1, UVISTA-Y-5, UVISTA-Y-6, UVISTA-J-1, and UVISTA-J-2) with plausible $z_{phot} \ge 8.5$ solutions, that stood out by their unprecedented brightness ($24.4 \le H \le 25.3$), which were detected in the UltraVISTA ultradeep stripes, and with coverage from the deepest optical ground-based data in that region to be followed up with HST/WFC3. Their positions and H-band fluxes are listed in Table 2.

### 3. HST Data and Photometry

The five bright candidate $z \sim 8$ sources presented in this work benefit from HST/WFC3 imaging obtained during the mid-cycle 24 (PI: R. Bouwens, PID: 14895). Observations

were performed from 2017 March 27 to 29. Table 2 summarizes the main observational parameters of the sample. Each source was observed for 1 orbit in total, subdivided as follows: $\sim 1500$ s ($\sim 0.65$ orbits) in the F098M band ($z_{098}$ hereafter), $\sim 460$ s (0.18 orbits) in the F125W band ($J_{125}$ hereafter), and $\sim 410$ s (0.17 orbits) in the F160W band ($H_{160}$ hereafter).

The field of UVISTA-Y-1 has also been observed by program 14808 (SUbaru Supernovae with Hubble Infrared—SUSHI; PI: Nao Suzuki) with $\sim 1000$ s integration time in the F105W band and $\sim 1200$ s integration time in the F140W band, which we included in our analysis. Image stamps in all the five HST bands centered at the position of UVISTA-Y-1 are shown in Figure 1.

The observations were processed using a customized version of multi-drizzle (Koekemoer et al. 2003). For each object, the images in the three HST bands F098M, F125W, and F160W were combined together into a red channel image. Figure 2 presents the image cutouts of the five objects in the three HST/WFC3 bands together with ground-based and Spitzer IRAC bands.

Photometry of the HST bands was extracted using SExtractor in dual image mode, with the detection performed on the red channel image. Fluxes were measured in apertures 0.6 wide (diameter) in each band, and corrected for the flux excluded by the finite aperture using the PSF curve of growth. The typical aperture corrections were $\le 3\%$ across the WFC3 bands, minimizing potential systematic effects from the different PSFs.

Using the new HST data, we also reprocessed the flux measurements in all the ground-based optical and NIR bands and in the Spitzer/IRAC bands. Fluxes were measured using the mophongo software in apertures 1.2 (diameter) and corrected to total using the light profile of each source. Remarkably, the optical data now include the mosaics from the HSC Survey (Aihara et al. 2017a, 2017b), not available at the time of the original selection of I. Labbé et al. (2017, in preparation). This new program provides deep observations in the g, r, i, z, and Y bands ($5\sigma$ depths of 26.6, 26.7, 26.2, 25.8, and 24.8 mag, respectively).

In Table 1 we list the average multiplicative corrections applied to convert aperture fluxes into total. For ground-based data they range from $\sim 1.6$ to $\sim 2.6$; for IRAC 3.6 $\mu m$ and 4.5 $\mu m$ data they are approximately $\sim 5.4$, while for the two reddest IRAC bands they have values of $\sim 8-10$. The NIR and IRAC bands are characterized by large aperture corrections,
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Image stamps, in inverted gray scale, centered at the position of UVISTA-Y-1, in the five *HST* bands available for this object, i.e., the three canonical bands targeted by our *HST* program (PI: R. Bouwens; *z*₀₉₈, *J*₁₂₅ and *H*₁₆₀) plus *Y*₁₀₅ and *JH*₁₄₀ from the SUSHI program (PI: N. Suzuki). The source is clearly detected in the *J*₁₂₅, *JH*₁₄₀ and *H*₁₆₀ bands, while it is only slightly detected (2.2σ) in *Y*₁₀₅. The cutouts in the *HST*/WFC3 bands have been smoothed with a σ = 1.4 pixel Gaussian kernel.</figcaption></figure>

which could introduce systematics in the estimates of the total
fluxes. As a sanity check, we repeated the photometry with an
aperture of 1"/8 diameter. The recovered total flux densities are
on average within a few percent of the measurements based on
the 1"/2 apertures, and within ~15% (~1σ) in the worst cases.
We therefore considered the photometry obtained adopting the
1"/2 aperture equally robust to that obtained with a larger
aperture, but with a higher S/N.

Uncertainties associated to flux densities were computed
differently depending on the data set. For HST/WFC3 bands,
we estimated the noise associated to the background from the
dispersion of values in 200 0.6 apertures randomly placed
across the image, free from sources according to the
segmentation map, and repeated this process 20 times to
increase the statistical significance. The final value was
obtained by applying the same aperture correction adopted
for the estimate of the total fluxes.

Uncertainties for ground-based and IRAC data were
computed by *mophongo*. Briefly, the rms of the pixels in
the residual image, obtained by subtracting all the detected
objects from the science cutout, was computed for apertures of
1.2. As an additional step, the rms value was taken to be the
maximum between the rms initially estimated by *mophongo*
and that obtained from the empty apertures (see Section 2 and
Table 1). The systematic errors of kernel reconstruction were
then added in quadrature and the result was scaled through the
aperture correction.

The uncertainties resulting from this method are therefore
not just the pure translation of the exposure map; specifically,
the introduction of the systematic error from the kernel
reconstruction and the scaling according to the aperture
correction, which itself is, in general, different for different
sources in a given band, makes the comparisons of uncertain-
ties across different sources in a given band less immediate.
Nonetheless, this method provides a robust and more
comprehensive estimate of flux uncertainties.

One example of the above behavior is given by the
uncertainties in the UltraVISTA bands of UVISTA-Y-5 and
UVISTA-Y-6. UVISTA-Y-6 lies at the border of one of the
ultradeep stripes, while UVISTA-Y-5 is located in the middle of
one of the ultradeep stripes. The ratio between the effective
exposure time of UVISTA-Y-6 to that of UVISTA-Y-5 is
~0.76, which would correspond to an increase of the rms
background for UVISTA-Y-6 by a factor of ~1.14. Instead, our
analysis recovers flux uncertainties higher for UVISTA-Y-5 than
for UVISTA-Y-6. Inspection of the *mophongo* output showed
that UVISTA-Y-5 is characterized by an rms background very

similar to that of UVISTA-Y-6 and by a larger aperture
correction (~2.5 versus ~2.0, suggesting a more extended
morphology for UVISTA-Y-5). The comparable values of the
rms background for the two sources are likely the result of a
larger value estimated by *mophongo* for the systematic
uncertainty associated to the kernel reconstruction for
UVISTA-Y-5.

As a further test, we repeated our analysis after replacing the
uncertainties with the maximum uncertainty measured for each
band across all the sources, and found results consistent with
those from the main analysis, further supporting our error
budget analysis.

The full set of measurements on ground- and space-based
mosaics are presented in Table 3. *HST* data were key to our
work as they provided accurate positional and morphological
priors for the *mophongo* photometry, enabling a more
accurate neighbor subtraction. This, together with the addi-
tional information provided by the HSC Survey (especially for
UVISTA-Y-1 which lacks coverage from the CFHTLS),
enabled a more accurate determination of the photometric
redshifts and stellar population parameters for the galaxies in
our sample.

4. Results: Improved Spectral Energy Distributions and
Photometric Redshifts

The left panels of Figures 3 and 4 present the spectral energy
distributions (SEDs) of the five sources studied in this work.
The filled green squares and arrows mark the WFC3
measurements and 2σ upper limits, respectively.

In order to further assess their nature, we computed
photometric redshifts running *EAZY* (Brammer et al. 2008)
with the standard set of SED templates together with three old-
and-dusty templates. Specifically, these templates correspond
to a 2.5 Gyr, single burst, passively evolving, $Z_\odot$, Bruzual &
Charlot (2003) stellar population, further reddened with
Calzetti et al. (2000) $A_V = 2.0, 5.0, 8.0$ mag curves.

The input catalog consisted of the flux measurements listed
in Tables 2 and 3. One of the advantages of working with flux
densities over magnitudes is that negative fluxes can be treated
in a natural way, without any need to convert them into upper
limits, thus preserving fidelity to observations.

Using the full set of bands, we find that UVISTA-Y-1,
UVISTA-Y-5, and UVISTA-Y-6 have photometric redshifts
$z_{phot} = 8.38_{-0.43}^{+0.35}$, $8.74_{-0.47}^{+0.45}$, and $8.53_{-0.80}^{+0.53}$, respectively with
$\chi^2 = 13.7, 10.6$, and 6.6. The remaining two sources
(UVISTA-J-1 and UVISTA-J-2) instead prefer solutions at
$z_{phot} \sim 2$ ($\chi^2 = 18.5$ and 23.2, respectively).

In Figure 3 we also present the best-fitting brown dwarf SED
template (light brown curve) and the best fit when we force the
solution to be at *z* < 6 (gray curve). Both these fits were
obtained considering the full set of photometric points. Neither
the brown dwarf nor the *z* < 6 solutions do a better job at
describing the observations compared to the *z* ~ 8 best-fit
template. Specifically, the brown dwarf template is inconsistent
with observations in the IRAC 3.6 and 4.5 µm bands and, for
UVISTA-Y-1, also with the *JH*₁₄₀ measurement. This is
reflected by the poorer best-fit *χ*², with *χ*² = 54.8, 58.5, and
41.2, respectively, for the three sources. Forcing the solution
to be at *z* < 6 results in *z*<sub>phot</sub> ~ 2. The best-fitting SED has a
substantial contribution from an old/dusty component and
provides a much better fit to the data than the brown dwarf
<table><caption>**Table 3**
Total Flux Densities for the Five Candidate $z \gtrsim 8$ LBGs over COSMOS/UltraVISTA Targeted by Our Small HST Program</caption><thead><tr><th>Filter</th><th>UVISTA-Y-1<br>(nJy)</th><th>UVISTA-Y-5<br>(nJy)</th><th>UVISTA-Y-6<br>(nJy)</th><th>UVISTA-J-1<br>(nJy)</th><th>UVISTA-J-2<br>(nJy)</th></tr></thead><tbody><tr><td>CFHTLS u*</td><td>...</td><td>-13 ± 17</td><td>8 ± 14</td><td>4 ± 15</td><td>13 ± 17</td></tr><tr><td>SSC B</td><td>-2 ± 7</td><td>4 ± 9</td><td>-11 ± 9</td><td>-4 ± 8</td><td>5 ± 11</td></tr><tr><td>HSC g</td><td>8 ± 15</td><td>-10 ± 20</td><td>4 ± 16</td><td>1 ± 15</td><td>6 ± 20</td></tr><tr><td>CFHTLS g</td><td>...</td><td>3 ± 15</td><td>1 ± 13</td><td>-8 ± 10</td><td>-4 ± 17</td></tr><tr><td>SSC V</td><td>-12 ± 17</td><td>-1 ± 24</td><td>0 ± 21</td><td>-7 ± 17</td><td>-11 ± 27</td></tr><tr><td>HSC r</td><td>-7 ± 14</td><td>-3 ± 18</td><td>11 ± 15</td><td>4 ± 15</td><td>10 ± 17</td></tr><tr><td>CFHTLS r</td><td>...</td><td>-15 ± 23</td><td>8 ± 21</td><td>-22 ± 17</td><td>7 ± 22</td></tr><tr><td>SSC r'</td><td>7 ± 16</td><td>-29 ± 23</td><td>20 ± 19</td><td>15 ± 16</td><td>-7 ± 18</td></tr><tr><td>CFHTLS y</td><td>...</td><td>-11 ± 28</td><td>22 ± 26</td><td>16 ± 23</td><td>-26 ± 28</td></tr><tr><td>CFHTLS i</td><td>...</td><td>-11 ± 29</td><td>5 ± 29</td><td>-28 ± 24</td><td>-29 ± 29</td></tr><tr><td>HSC i</td><td>21 ± 21</td><td>-20 ± 27</td><td>1 ± 23</td><td>19 ± 21</td><td>23 ± 27</td></tr><tr><td>SSC i<sup>+</sup></td><td>-9 ± 22</td><td>-36 ± 26</td><td>2 ± 21</td><td>-18 ± 24</td><td>-24 ± 27</td></tr><tr><td>CFHTLS z</td><td>...</td><td>6 ± 63</td><td>-13 ± 59</td><td>-65 ± 52</td><td>32 ± 62</td></tr><tr><td>HSC z</td><td>9 ± 31</td><td>-27 ± 39</td><td>17 ± 33</td><td>52 ± 31</td><td>8 ± 37</td></tr><tr><td>SSC z'</td><td>51 ± 64</td><td>-51 ± 93</td><td>69 ± 85</td><td>70 ± 64</td><td>-39 ± 83</td></tr><tr><td>HSC Y</td><td>-31 ± 76</td><td>-91 ± 98</td><td>89 ± 80</td><td>80 ± 76</td><td>-29 ± 98</td></tr><tr><td>z<sub>098</sub></td><td>46 ± 34</td><td>-7 ± 34</td><td>9 ± 34</td><td>13 ± 34</td><td>22 ± 34</td></tr><tr><td>UVISTA Y</td><td>18 ± 48</td><td>-42 ± 68</td><td>16 ± 51</td><td>67 ± 55</td><td>37 ± 67</td></tr><tr><td>Y<sub>105</sub></td><td>92 ± 41</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>J<sub>125</sub></td><td>279 ± 70</td><td>195 ± 70</td><td>172 ± 70</td><td>220 ± 70</td><td>304 ± 70</td></tr><tr><td>UVISTA J</td><td>324 ± 50</td><td>235 ± 66</td><td>211 ± 53</td><td>125 ± 50</td><td>195 ± 66</td></tr><tr><td>JH<sub>140</sub></td><td>303 ± 44</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>H<sub>160</sub></td><td>511 ± 107</td><td>152 ± 107</td><td>272 ± 107</td><td>188 ± 107</td><td>257 ± 107</td></tr><tr><td>UVISTA H</td><td>455 ± 61</td><td>393 ± 86</td><td>280 ± 66</td><td>510 ± 70</td><td>657 ± 86</td></tr><tr><td>UVISTA K<sub>S</sub></td><td>480 ± 77</td><td>321 ± 102</td><td>271 ± 82</td><td>602 ± 95</td><td>822 ± 102</td></tr><tr><td>IRAC 3.6 µm</td><td>623 ± 85</td><td>289 ± 74</td><td>434 ± 106</td><td>1162 ± 90</td><td>913 ± 110</td></tr><tr><td>IRAC 4.5 µm</td><td>931 ± 109</td><td>589 ± 86</td><td>598 ± 130</td><td>1277 ± 102</td><td>1204 ± 130</td></tr><tr><td>IRAC 5.8 µm</td><td>-2893 ± 2568</td><td>-1978 ± 4831</td><td>-643 ± 3000</td><td>-2209 ± 2655</td><td>8075 ± 4602</td></tr><tr><td>IRAC 8.0 µm</td><td>1423 ± 3021</td><td>499 ± 6310</td><td>-3325 ± 3803</td><td>1992 ± 3771</td><td>6734 ± 5613</td></tr></tbody></table>

Note. Measurements for the ground-based and Spitzer/IRAC bands are 1′′2 aperture fluxes from mophongo corrected to total using the PSF and luminosity profile information; HST/WFC3-band flux densities are measured in 0′′6 apertures and converted to total using the PSF growth curves. Flux density measurements of three other ultra-luminous ($M_{UV} \lesssim -22$) $z \sim 8$ candidates are presented in Table 6 from Appendix B.

However, it remains in tension with the data, resulting in $\chi^2$ values of $\chi^2 = 17.6, 25.3$, and $8.7$, respectively, for the three sources, i.e., $\Delta\chi^2 = 3.9, 14.7$, and $2.1$, respectively, worse than the $z \sim 8$ fits. Similarly, the best-fit brown dwarf templates for UVISTA-J-1 and UVISTA-J-2, displayed in Figure 4, are inconsistent with our measured fluxes in the IRAC 3.6µm and 4.5µm bands, where $\chi^2 = 61.3$ and 71.9, respectively.

To ensure that our photometric redshift results are robust against potential errors and underestimates of the flux uncertainties for individual sources, we perturbed these by factors 1–1.5 randomly extracted from a uniform distribution. The new catalog was analyzed following our standard procedure and the whole process was repeated 500 times. All of the recovered best-fit redshifts were within the $1\sigma$ uncertainties of our nominal $z \sim 8$ solutions.

We also looked at what happened to our photometric redshift solutions if the flux uncertainties were somewhat smaller than what we estimate, as for example we found in Section 2 (typical differences were factors of 1.5). We found photometric redshifts $z_{phot} = 8.13^{+0.38}_{-0.42}, 8.57^{+0.47}_{-0.47}$, and $8.43^{+0.57}_{-0.89}$, for UVISTA-Y-1, UVISTA-Y-5, and UVISTA-Y6, respectively, with associated probabilities for the solution to be at $z > 6$ of $p(z > 6) = 0.99, 0.99, \text{and } 0.84$, respectively.

Limitations in our knowledge of the intrinsic SED shapes of $z \gtrsim 7$ galaxies (e.g., Balmer break amplitude, nebular emission

lines equivalent width) make fits to the redder wavelength data more difficult, particularly in our attempts to derive accurate redshifts for the sources. During the SED fitting process, non-null colors from contiguous broad bands can be misinterpreted as features which are not intrinsic to the source under analysis. For this reason, we repeated the photometric redshift measurements of the three $z \gtrsim 8$ sources after excluding the IRAC 3.6, 4.5, and 5.8 µm bands, as these are likely contaminated by strong emission lines and/or potentially contain the Balmer/ 4000 Å break. Both these properties are still poorly determined at these redshifts and any assumption about them could therefore introduce systematics in the redshift estimates. However, we still included the 8.0 µm data as they are likely not contaminated by strong nebular emission, yet provide constraints for the SED modeling. Indeed, the $K_S - [4.5]$ color could be interpreted as the Balmer break, guiding the fit toward higher-redshift solutions. These new measurements resulted in lower photometric redshifts: $z_{phot} = 8.02^{+0.41}_{-0.49}, 8.39^{+0.60}_{-0.60}$, and $8.35^{+0.65}_{-0.81}$ for UVISTA-Y-1, UVISTA-Y-5, and UVISTA-Y-6, respectively. For this reason, we consider our fiducial photometric redshifts for UVISTA-Y-1, UVISTA-Y-5, and UVISTA-Y-6 those obtained without the IRAC bands. We remark here, however, that the IRAC bands are nevertheless useful for our interpretation of these sources as they allow us to distinguish between genuine high-redshift sources and lower-redshift interlopers.
The photometric redshift measurements for UVISTA-J-1 and UVISTA-J-2 were repeated after excluding the HST/WFC3 and HSC flux measurements, to explore the possible reasons for the detected change in redshift. The redshift of UVISTA-J-2 obtained without the WFC3 and HSC bands is $z_{phot} = 10.1_{-0.8}^{+1.4}$, consistent with the initial selection. The new $H_{160}$ observations point to a much redder overall NIR color (e.g., $H_{160}$ - [3.6]) for the source, indicating that the $z < 6$ solution is clearly the best one. For UVISTA-J-1, however, the photometric redshift we find does not sensibly change ($z_{phot} = 2.2_{-0.4}^{+0.6}$). After further inspection, we conclude that the likely reason for this is a higher flux measurement in the 3.6 µm band we obtained using the new HST data set as morphological and positional prior which allowed for a more accurate subtraction of the neighbors, compared to the initial estimate obtained adopting the UltraVISTA bands. Both cases further stress the importance of high-resolution imaging from HST in ascertaining the nature of candidate high-z sources.

None of the five sources has a counterpart in the deep VLA catalogs of Smolcic et al. (2017) nor in the X-ray catalogs from XMM and Chandra (Cappelluti et al. 2009; Civano et al. 2016; Marchesi et al. 2016, respectively). Finally, visual inspection of the MIPS µm mosaic from the S-COSMOS project (Sanders et al. 2007) did not show any evidence for the presence of sources at the nominal positions; we note, however, that a source is likely present ~0.8 east of UVISTA-J-1.

The best-fit SEDs are shown as solid curves in the plot of the left-side panels of Figures 3 and 4, while the right-side panels show the redshift likelihood generated by EAzY. In the following paragraphs we comment on the individual sources.

**UVISTA-Y-1:** This source is undetected (<2σ) in the HST/WFC3 $z_{098}$ band, strengthening the evidence that this is a $z > 7.5$ LBG. The WFC3 photometry in the $J_{125}$ and $H_{160}$ is consistent at 1σ or better with that in the UltraVISTA J and H bands, respectively. This source also benefits from additional WFC3 coverage in the $Y_{105}$ and $JH_{140}$ bands from the SUSHI program (PI: N. Suzuki). The measurement in the $JH_{140}$ band is consistent with the best-fit SED. The $p(z)$ is characterized by a solution at $8.02_{-0.49}^{+0.41}$, with a marginal secondary peak at $z \sim 1.8$ ($p(z < 6) = 0.12$). The $Y_{105}$ band shows a 2.2σ detection, as expected if the source is at $z \sim 8$.

**UVISTA-Y-5:** The source is undetected (<2σ) in the HST/WFC3 $z_{098}$ band, strongly favoring a $z > 7.5$ solution for this source. The WFC3 photometry in the $J_{125}$ is consistent at 1σ with that in the UltraVISTA J band. However, the flux measurement in the $H_{160}$ results in a 1.42σ detection only, and is consistent with the UltraVISTA H-band photometry at ~3σ level. In Appendix A we analyze more in detail the main effects that could explain the systematic differences observed in the $H_{160}$ and UltraVISTA H. Here we caution the reader that the observed discrepancy reduces our confidence on the high-redshift solution. The $p(z)$ is characterized by a peaked distribution at $8.39_{-0.60}^{+0.60}$, with a very marginal secondary peak at $z \sim 1.8$ ($p(z < 6) = 0.009$).

**UVISTA-Y-6:** This source is undetected (<2σ) in the HST/WFC3 $z_{098}$ band, again favoring a $z > 7.5$ solution. The WFC3 photometry in the $J_{125}$ and $H_{160}$ is consistent at 1σ or better with that in the UltraVISTA J and H bands, respectively. The $p(z)$ shows a distribution with best-fit solution $z_{phot} = 8.35_{-0.81}^{+0.65}$ with a hint of secondary solution at $z < 6$ ($p(z < 6) = 0.168$).

**UVISTA-J-1:** The object is formally undetected in the $H_{160}$ band (<2σ), making it consistent with the UltraVISTA H band

only at ~3.5σ level. We again refer the reader to Appendix A for a detailed discussion about possible origins of the observed difference, and flag this source because of the decreased confidence on the redshift determination. The photometry in the $J_{125}$ band, instead, is consistent with that in the UltraVISTA J band at 1σ level. The lower-z solution is enforced by the fact that the source is detected in the HSC z band at 1.7σ level. The fiducial photometric redshift is $2.05_{-0.46}^{+0.49}$. The $p(z)$ shows a peaked distribution around $z \sim 2$ with no further secondary peaks at higher redshifts.

**UVISTA-J-2:** The HST WFC3 photometry in the $J_{125}$ and $H_{160}$ is consistent with that in the corresponding UltraVISTA bands at ~1.8–3σ, respectively, with nominal redshift of $1.96_{-0.33}^{+7.04}$. Similarly to what was done for UVISTA-Y-5 and UVISTA-J-1, in Appendix A we analyze the main effects that could generate the observed difference between the $H_{160}$ and H bands. Again, here we caution the reader that this discrepancy reduces our confidence on our redshift estimate. The large uncertainty associated to the upper limit makes it consistent with $z \sim 9–9.5$; however, the $p(z)$ shows a pronounced peak at $z \sim 2$ and a secondary peak at $z \sim 10$, with a likelihood for the SED to be at $z > 6$ of $p(z > 6) = 0.066$. For this reason we consider the fiducial redshift for this source to be the $z \sim 2$ solution.

# 5. Discussion

## 5.1. The Brightest Candidate LBGs at $z \gtrsim 8$

Figure 5 presents our sample of candidate $z \gtrsim 8$ LBGs in the redshift-$M_{UV}$ plane, together with recent LBG selections covering the bright end of the UV LF at $z \gtrsim 7$ of Oesch et al. (2014) on XDF/HUDF, Bouwens et al. (2015), Roberts-Borsani et al. (2016), and Oesch et al. (2016) based on CANDELS data, Calvi et al. (2016) from the BoRG program, Bowler et al. (2017) from UltraVISTA, and Ono et al. (2017) from the HSC survey. We note, however, that the candidates of Calvi et al. (2016) lack IRAC coverage and those of Ono et al. (2017) have measurements only at optical wavelengths from the HSC Survey, resulting in their nature being more uncertain. The three $z \sim 8$ galaxies reported on here constitute the brightest, most reliable $z \sim 8–9$ galaxies discovered to date. In order to put their luminosities in better context, in the same figure we also represent the evolutionary relation of the characteristic magnitude of the UV LF of Bouwens et al. (2015) up to $z = 8$ and its extrapolation to $z \sim 10$. Our sample of luminous galaxies are ~1.8 mag more luminous (a factor ~5.3×) than the estimated characteristic magnitude at $z \sim 9$.

## 5.2. $\beta - M_{UV}$ Relation

We measured the rest-frame UV slope ($\beta$) by fitting a power law of the form $f_{\lambda} \propto \lambda^{\beta}$ to the fluxes in the $H_{160}$ and in the UltraVISTA H and $K_S$ bands. The results are presented in Figure 6 and listed in Table 4. These slopes have an inverse-variance-weighted average value of $\beta = -1.74 \pm 0.35$ and are consistent with the recent determination of the UV slopes of Bowler et al. (2017) for LBGs with similar luminosity ($M_{UV} \sim -22.5$) identified at $z \sim 7$ over the COSMOS/ UltraVISTA field, suggesting a slow evolution of $\beta$ for luminous galaxies at early cosmic epochs. Our measurements are also consistent with the UV slope $\beta = -2.1 \pm 0.3$ from stacking of bright ($M_{UV} \sim -21$) $z \sim 10$ LBGs by Wilkins et al. (2016). For comparison, in the plot we also show the bi-
<table><caption>**Table 4**

Physical Parameters for the Three Galaxies with Photometric Redshift $z \gtrsim 8$<sup>a</sup></caption><thead><tr><th>ID</th><th>z<sub>photo</sub><sup>b</sup></th><th>p(z > 6)<sup>c</sup></th><th>M<sub>UV</sub><sup>d</sup><br>(mag)</th><th>β<sup>e</sup></th><th>EW<sub>0</sub>(Hβ + [O III])<sup>f</sup><br>(Å)</th></tr></thead><tbody><tr><td>UVISTA-Y-1</td><td>8.02<sup>+0.41</sup><sub>-0.49</sub></td><td>0.88</td><td>-22.46 ± 0.15</td><td>-1.98 ± 0.67</td><td>1041<sup>+713</sup><sub>-515</sub></td></tr><tr><td>UVISTA-Y-5<sup>g</sup></td><td>8.39<sup>+0.60</sup><sub>-0.60</sub></td><td>0.99</td><td>-22.37 ± 0.24</td><td>-1.62 ± 0.42</td><td>887<sup>+1323</sup><sub>-686</sub></td></tr><tr><td>UVISTA-Y-6</td><td>8.35<sup>+0.65</sup><sub>-0.81</sub></td><td>0.83</td><td>-21.97 ± 0.26</td><td>-2.12 ± 1.44</td><td>1291<sup>+1749</sup><sub>-940</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> The main properties of three other ultra-luminous ($M_{UV} \lesssim -22$) $z \sim 8$ candidates are presented in Table 7 from Appendix B.

<sup>b</sup> Best photometric redshift estimate from EAzY, excluding the IRAC bands from the fit, and corresponding 68% confidence interval.

<sup>c</sup> Probability, computed by EAzY, that the solution is at $z > 6$.

<sup>d</sup> Absolute magnitudes at rest-frame 1600 Å from EAzY.

<sup>e</sup> Rest-frame UV continuum slopes from the H<sub>160</sub> and UltraVISTA H and K<sub>S</sub> bands.

<sup>f</sup> Rest-frame equivalent width of Hβ + [O III] obtained from the K<sub>S</sub> - [4.5] color assuming an SED flat in $f_ν$ (i.e., $β = -2$). If $β = -1.8$, the EW and associated uncertainties would be a factor ~1.2 smaller.

<sup>g</sup> Given tension in the flux measurement between the H<sub>160</sub> and UltraVISTA H bands, we caution the reader that the redshift estimate for this object (and therefore M<sub>UV</sub>, β, and EW<sub>0</sub>) is less robust than that of the other two sources in this table. See Section 4 and Appendix A for further details.

weight UV slope measurements at $z \sim 7$ and $z \sim 8$ from Bouwens et al. (2014), using data from the CANDELS GOODS-N, CANDELS GOODS-S, and the HST HUDF/XDF fields.

Recent works have identified a correlation between the UV luminosity and the slope of the UV continuum and as a function of redshift: redder slopes are observed at fixed redshift for more luminous galaxies and at fixed luminosity for galaxies at later cosmic times (Wilkins et al. 2011; Bouwens et al. 2012, 2014; Castellano et al. 2012; Finkelstein et al. 2012; Dunlop et al. 2013; Jiang et al. 2013; Duncan et al. 2014; Rogers et al. 2014; Duncan & Conselice 2015; see also Oesch et al. 2013 and Stefanon et al. 2017 for similar relations of β and rest-frame optical luminosities). This behavior has been interpreted as the emergence of older stellar populations, dust, and metals in more luminous galaxies. In Figure 6 we also plot the recent determination of the β-M<sub>UV</sub> relation of Bouwens et al. (2014) at $z \sim 7$ and $z \sim 8$. Our measurements lie below their extrapolation to the luminosity range probed in this work, although the large uncertainties associated to our β measurements make them consistent at $<1σ$ with those relations, thus preventing us from further inspecting any differential evolutionary path of β with luminosity and redshift.

## 5.3. Size Measurement

The availability of high-resolution imaging from our small HST/WFC3 program allowed us to pursue a first study of the size and morphological properties of extremely bright $z \sim 8$ galaxies.

Morphological information was recovered by running galfit (Peng et al. 2002, 2010), which fits the convolution of a brightness profile with a PSF. The advantage of this approach is that the extracted morphological parameters are deconvolved from PSF effects. For this work, we considered only the Sèrsic (1968) profile, characterized by an effective radius $R_e$ and an index (n) expressing how steeply the wings of the profile decrease with the radius. We note here that the symmetry of the brightness profile assumed by the Sèrsic form could result in an over-simplification, and consequently limitations, at the time of describing the morphological properties of high-redshift galaxies (e.g., $R_e$ in the presence of clumpy or merging systems). Indeed, recent studies have shown that sources are observed to be non-symmetric over a

wide range of redshifts (e.g., Law et al. 2007; Mortlock et al. 2013; Huertas-Company et al. 2015; Ribeiro et al. 2016; Bowler et al. 2017), suggesting that high-redshift galaxies could be characterized by a range of sizes and morphologies, resulting from different physical processes.

Considering that the limited S/N of our observations does not allow us to perform a more comprehensive and detailed morphological analysis, we base our analysis on the working hypothesis of a symmetric Sèrsic profile. Furthermore, because of the relatively low S/N in most of the WFC3 images, for our analysis we only considered the JH<sub>140</sub> band of UVISTA-Y-1 (~rest-frame 1600 Å), i.e., the highest S/N observation for the brightest object.

The first estimate of the target position, its magnitude, its $R_e$, the axis ratio, and the value of the local background, needed as input by galfit, was obtained from SExtractor. During the fitting process, we left the $R_e$, magnitude, and axis ratio free to vary, while we kept the background fixed. Because of the small extension of the brightness profile and of the relatively low S/N of our data, during the fitting process we fixed the Sèrsic index to $n = 1.5$, consistent with measurements at $z \sim 7-10$ (e.g., Oesch et al. 2010; Holwerda et al. 2015; Bowler et al. 2017). We then verified that $R_e$ does not systematically change (≤10%) when the Sèrsic index varies in the range $1.2 < n < 2.0$. This variation was added in quadrature to the uncertainty on $R_e$ provided by galfit. In order to ensure the most robust result, in the fit we also included all the neighbors within 5"/0 from the nominal position of UVISTA-Y-1. Because the $R_e$ directly provided by galfit corresponds to the major semi-axis, and in order to compare to estimates from the literature, we circularized it as $R_{e,circ} = R_e \sqrt{b/a}$, where $b/a$ is the minor-to-major axis ratio. As a consistency check, we also derived $R_e$ using SExtractor. In this case, the final value for $R_{e,SE} = \sqrt{R_{e,obs}^2 + r_{PSF}^2}$, with $R_{e,obs}$, SE the effective radius measured by SExtractor and $R_{PSF}$ that of the JH<sub>140</sub> PSF, with $R_{PSF} = 0.12$.

We find $R_{e,circ} = 0.9 \pm 0.3$ kpc from galfit, consistent with $R_{e,SE} = 0.7$ kpc estimated with SExtractor. In Table 5 we list the main morphological parameters we obtain from the two methods. Our values are consistent at $1σ$ with estimates of $R_e$ for $M_{UV} \sim -22$ LBGs at $z \sim 7$ ($R_{e,circ} = 0.6 - 0.9 \pm 0.2$ from a stacking analysis; Bowler et al. 2017) and $z \sim 11$.
<table><caption>**Table 5**

Morphological Parameters for UVISTA-Y-1 Measured on the JH<sub>140</sub> Band with galfit and SExtractor</caption><thead><tr><th>Algorithm</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>R<sub>e,circ</sub><sup>a</sup><br>(kpc)</th><th>q<sup>b</sup></th><th>n<sup>c</sup></th><th>Σ<sub>SFR</sub><sup>d</sup><br>(M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup>)</th></tr></thead><tbody><tr><td>galfit</td><td>09:57:47.910</td><td>+02:20:43.50</td><td>0.9 ± 0.3</td><td>0.9 ± 0.2</td><td>1.5</td><td>11<sup>+14</sup><sub>-5</sub></td></tr><tr><td>SExtractor</td><td>09:57:47.910</td><td>+02:20:43.50</td><td>0.7 ± 0.1</td><td>0.6 ± 0.1</td><td>...</td><td>17<sup>+10</sup><sub>-5</sub></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Circularized effective radius.

<sup>b</sup> Minor-to-major axis ratio.

<sup>c</sup> Sérsic index. This was kept fixed when running galfit.

<sup>d</sup> Star formation rate surface density, computed following Ono et al. (2013).

$(R_{e,circ} = 0.6 \pm 0.3)$ for the brightest known galaxy at the highest redshift, with luminosity similar to that of our sample; Oesch et al. 2016). Moreover, because the evolution of the characteristic luminosity of the UV LF is small for $z \sim 4-10$ (Bouwens et al. 2015; Finkelstein et al. 2015; Bowler et al. 2017; Ono et al. 2017), and considering that our sources constitute the very bright end of the UV LF, the absolute magnitude corresponding to a constant cumulative number density should evolve very little over $z \sim 4-10$. This means that, under the further assumption of a smooth evolution of the star formation history (SFH), selecting galaxies with approximately the same (high) luminosity corresponds to selecting the descendants of the luminous galaxies observed at the highest redshift in the sample.

In the top panel of Figure 7 we present a compilation of size measurements for LBGs at $z > 4$ and $M_{UV} \sim -22$ from Huang et al. (2013), Shibuya et al. (2015), Oesch et al. (2016), Curtis-Lake et al. (2016), and Bowler et al. (2017). The plot suggests only a modest evolution in size for luminous galaxies (factor of $\sim 3 \times$) during approximately the first 1.5 Gyr of cosmic time. The bottom panel of Figure 7 presents the evolution of the SFR surface density ($\Sigma_{SFR}$), computed using the recipe of Ono et al. (2013). The SFR is estimated from the UV luminosity following Kennicutt (1998) under the assumption of negligible dust obscuration. The SFR is then divided by the area corresponding to $R_{e,circ}$ and applying a further factor 0.5 to take into account that observationally we can only access approximately half of the surface of each galaxy. The value we find for $\Sigma_{SFR} \sim 11 M_{☉} \text{yr}^{-1} \text{kpc}^{-2}$ is consistent with measurements at lower luminosities (e.g., Ono et al. 2013; Holwerda et al. 2015; Shibuya et al. 2015, although Oesch et al. 2010 found $\Sigma_{SFR}$ for $L \lesssim L^*$ galaxies a factor $\sim 3 \times$ lower).

Interestingly, but unsurprisingly, $\Sigma_{SFR}$ decreases with cosmic time, although with marginal statistical significance. Some recent studies of $z \sim 4-8$ LBGs have found indication for a non-evolving $\Sigma_{SFR}-z$ relation (e.g., Oesch et al. 2010; Ono et al. 2013). This is qualitatively consistent with the increase of the SFR density with cosmic time combined with the increase in size. Our (mildly) evolving $\Sigma_{SFR}$, instead, is the direct consequence of the evolution in size of galaxies with luminosity approximately constant over $4 \lesssim z \lesssim 10$.

We finally note that recent methods for the morphological analysis of high-redshift galaxies have found that the evolution of size could have been much less pronounced than recovered through more classical approaches (e.g., Law et al. 2007; Curtis-Lake et al. 2016; Ribeiro et al. 2016). While there is no reason to exclude that this could be the case at even higher redshift, data with better S/N are necessary for a more robust assessment.

**5.4. Volume Density at $z \sim 8$**

Using the results obtained in the previous sections, we estimate the contribution of the three candidate $z \sim 8$ LBGs to the UV LF. Here we focus on the HST sample analyzed in this work, which constitute the brightest end of the UV LF. A more comprehensive UV LF including the complete sample of fainter sources detected over COSMOS/UltraVISTA will be presented in I. Labbé et al. (2017, in preparation). The measurement of the volume density relies on estimating the detection completeness and the selection function associated to our selection criteria. We recovered these two quantities using similar procedures as described in Bouwens et al. (2015). Briefly, 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 re-sized to account for the change in angular diameter distance with redshift and for evolution of galaxy sizes at fixed luminosity (effective radius $r_e \propto (1+z)^{-1}$: Oesch et al. 2013; Ono et al. 2013; Holwerda et al. 2015; Shibuya et al. 2015). The template images were 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 and consistent with the determinations from this work (Bouwens et al. 2012; Finkelstein et al. 2012; Rogers et al. 2014). The simulations included 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 assumed a continuum flat in $f_ν$ and emission lines with fixed rest-frame EW(Hα+[N II]+[S II]) = 300 Å and rest-frame EW([O III] + Hβ) = 500 Å, consistent with the results of Labbé et al. (2013), Stark et al. (2013), Smit et al. (2014, 2015), and Rasappu et al. (2016). The same detection and selection criteria as described in Section 2 were then applied to the simulated images to recover the completeness as a function of magnitude and redshift.

Given that the source detection was performed on the UltraVISTA mosaics, roughly characterized by a dual depth (ultradeep and deep), the above process was independently executed in regions corresponding to the two depths. Figure 8 presents the selection functions associated to our criteria for the UltraVISTA ultradeep and deep stripes, used to estimate the comoving volumes entering the LF determinations. The plots show that in the ultradeep stripes our criteria allow us to select galaxies at $7.1 \lesssim z \lesssim 9.2$. In the deep stripes, instead, the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Image stamps in inverted gray scale of the five bright candidate $z \gtrsim 8$ LBGs in the stacked optical, stacked Y, HST/WFC3 $z_{098}$, $J_{125}$, $H_{160}$, UltraVISTA J, H, and $K_s$ and Spitzer IRAC 3.6 and 4.5µm bands. Each row corresponds to a source, as labeled on the left, where we omitted the prefix UVISTA from the object name for clarity. Each cutout is 5′′0 × 5′′0. The ground-based and IRAC cutouts are shown after removing the neighbors. The cutouts in the HST/WFC3 bands have been smoothed with a $\sigma = 1.4$ pixel Gaussian kernel. Postage stamp images of three other ultra-luminous ($M_{UV} \lesssim -22$) $z \sim 8$ candidates are presented in Figure 10 from Appendix B.</figcaption></figure>

range of redshift selection is slightly broader, $6.9 \lesssim z \lesssim 9.3$, qualitatively consistent with the fact that shallower depths in the NIR bands can also accommodate slightly different solutions.

The volume density associated to the three $z \gtrsim 8$ candidate LBGs was computed using the $1/V_{max}$ method (Schmidt 1968), and following the prescription of Avni & Bahcall (1980) for a coherent analysis, in order to deal with the different depths of the deep and ultradeep stripes of the UltraVISTA field. The $1/V_{max}$ method is intrinsically sensitive to local overdensities of galaxies; however, given the small sample considered in this work, we consider its potential effects by including the cosmic variance in the error budget. On the other hand, the $1/V_{max}$ method directly provides the normalization of the LF.

Considering that the absolute magnitudes of the three $z \gtrsim 8$ candidate LBGs are within 0.5 mag, the volume density was computed in one bin only. We obtained a volume density of $\Phi = 8.49_{-4.60}^{+8.23} \times 10^{-7} \text{ Mpc}^{-3} \text{ mag}^{-1}$ at $M_{UV} = -22.21 \pm 0.25$. The uncertainties associated to the volume density were computed following the recipe of Gehrels (1986), and adding in quadrature 24% of cosmic variance following Moster et al. (2011). Our measurement is shown in the top panel of Figure 9 with a filled red circle, together with a compilation of previous determinations of the bright end of the UV LF at $z \sim 8$. To avoid potential systematics, we limit our comparison to studies based on field galaxies (Bradley et al. 2012; McLure et al. 2013; Schenker et al. 2013; Schmidt et al. 2014; Bouwens et al. 2015; Finkelstein et al. 2015), and exclude UV LFs from studies based on galaxy cluster fields. Our measurement constitutes the first volume density estimate for $M_{UV} \lesssim -22$ at $z \sim 8$ with confidence $\gtrsim 1\sigma$ and is consistent with previous upper limits. In the same panel we also reproduce the Schechter

(1976) parameterization of the UV LF at $z \sim 8$ from Bouwens et al. (2015). Our estimate of the bright end agrees well with the exponential decline of the current Schechter form. Since the sample of I. Labbé et al. (2017, in preparation) includes three more potential galaxies at $7.5 \lesssim z_{phot} \lesssim 8.5$, which however did not enter the selection criteria for the HST proposal, here we also present the volume density obtained including all the six sources. The multi-wavelength photometry and results of SED fitting for these three additional sources are presented in Appendix B. At $M_{UV} = -22.21 \pm 0.25$ the volume density is $\Phi = 17.3_{-6.9}^{+10.3} \times 10^{-7} \text{ Mpc}^{-3} \text{ mag}^{-1}$. This measurement is plotted in Figure 9 with a pink filled circle and it is still consistent with the recent determinations of the bright end of the $z \sim 8$ UV LF(s).

Recent studies of the bright end of the UV LF at $z \gtrsim 6$ suggest that the LF could be parameterized by a double power law (DPL; Bowler et al. 2014, 2015, 2017; Ono et al. 2017) originated by an excess of luminous galaxies compared to the exponential decline of the Schechter function. The magenta curve in Figure 9 presents the DPL of Bowler et al. (2017) after evolving the faint-end slope, characteristic magnitude, and normalization factor to $z \sim 8$ using the evolution of Bouwens et al. (2015, see their Section 5.1). The DPL well describes the points, in particular considering the measurements of McLure et al. (2013) at absolute magnitudes brighter than $L^*$. Our single measurement is not able to distinguish between the two scenarios, however, as the corresponding absolute magnitude lies at the intersection of the Schechter and the DPL forms.

Because the $z_{phot}$ solutions for half of our sample of candidate $z \gtrsim 8$ LBGs may have values close to $z_{phot} \sim 9$ when including the IRAC bands, we also computed the volume density associated to the three sources (UVISTA-Y-5,
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left panels: SED for the three LBGs at $z \gtrsim 8$. The colored squares with black errorbars mark the photometric measurements, while arrows represent $2\sigma$ upper limits. Open squares and arrows mark the IRAC 3.6, 4.5, and 5.8µm bands, not used for the measurement of the fiducial photometric redshift. Photometry in the HST/WFC3 bands is indicated by the green points and arrows; HSC Survey data is represented in yellow. The fiducial best fit SED template from EAzY is indicated by the thick blue curve; the thin dark blue curve represents the best-fit SED when all bands are used for the photometric redshift measurement; the light brown curve presents the best-fitting brown dwarf template, while the gray curve the solution obtained forcing the redshift to be $z < 6$. Right panels: redshift likelihood distributions ($p(z)$) for the three LBGs for the fiducial solution (blue) and for the solution obtained considering the full set of flux measurements. The label in the top-left corner indicates the estimated photometric redshifts. The $p(z)$ are peaked, with no or very low integrated probability for a secondary solution at lower redshifts. We caution the reader, however, that given the flux inconsistency between the $H_{160}$ and the UltraVISTA *H* bands, the redshift estimate for UVISTA-Y-5 may be less robust than that of the other two sources; further details are discussed in Section 4 and Appendix A. The SEDs of three other ultra-luminous ($M_{UV} \lesssim -22$) $z \sim 8$ candidates are presented in Figure 11 from Appendix B.</figcaption></figure>

UVISTA-Y-6, and UVISTA-Y-2) with $z_{phot} > 8.5$. We obtain
$\Phi = 8.5_{-4.6}^{+8.2} \times 10^{-7} \text{ Mpc}^{-3} \text{ mag}^{-1}$. In the lower panel of
Figure 9 we compare our $z \sim 9$ volume density measurement
with the UV LFs at $z \sim 9$ from Oesch et al. (2013), McLure
et al. (2013), Bouwens et al. (2016), Calvi et al. (2016),
McLeod et al. (2016), and Ishigaki et al. (2017). Our estimate is
consistent with the measurement of Calvi et al. (2016),
although it corresponds to higher densities than expected from
the Schechter determination of Bouwens et al. (2016). In the
same panel we also plot the bright end of the DPL we
constructed for the $z \sim 8$ bin, renormalized to match the
density of the Schechter form at the characteristic magnitude. It
agrees within the error bars with our volume measurement.

The volume density we estimate at $z \sim 8$ is consistent with
that at $z \sim 9$, albeit with large statistical uncertainties, and
suggests a slow evolution of the brightest objects at early

cosmic epochs. Remarkably, this is still valid considering that
our volume density measurements are consistent at $\sim 1\sigma$ with
that for the $z \sim 11$ source GN-z11 (Oesch et al. 2016).
Assuming a smooth SFH, this could imply that these bright
(and possibly massive) galaxies assembled extremely rapidly in
the first few hundred Myr after the big bang. A very bursty
SFH, instead, would make any interpretation challenging,
because the number density would be a (random) combination
of bright (massive) galaxies with reduced SFR and lower-mass
galaxies with strong SFRs.

**6. Summary and Conclusions**

Here we report on HST/WFC3/IR observations on five very
bright $z \sim 8-9$ candidates identified over UltraVISTA. The
targeted sources were drawn from a sample of 16 very bright
$z \sim 8-9$ galaxies identified by I. Labbé et al. (2017, in
<figure><img src="image_7.png" /><figcaption>**Figure 4.** SEDs and redshift likelihood distributions for the two galaxies with likely best solution at $z \sim 2$. The result of excluding HST/WFC3 and HSC measurements for the recovery of the photometric redshifts is marked by the gray curves. Other plotting conventions as in Figure 3. We caution the reader, however, that given the flux inconsistency between the $H_{160}$ and the UltraVISTA H bands, the redshift estimates for these two sources may be less robust; further details are discussed in Section 4 and Appendix A.</figcaption></figure>

preparation), and constituted the brightest sources from that sample ($24.5 < H < 25.2$) with a plausible $z \gtrsim 8.5$ solution. The five sources in this sample (labeled UVISTA-Y-1, UVISTA-Y-5, UVISTA-Y-6, UVISTA-J-1 and UVISTA-J-2) that stood out for their brightness ($24.5 \lesssim H \lesssim 25.2$), for having plausible redshift $z_{phot} \sim 8.5$ solutions, for being positioned over the UltraVISTA deep stripes and had coverage from the deepest optical ground-based data, have recently been observed with HST/WFC3 (PI: R. Bouwens, PID: 14895) to try to confirm their nature.

The present work is devoted to the analysis of those sources specifically targeted with HST/WFC3 follow-up observations. Nevertheless, this analysis does present three other ultra-luminous $z \sim 8$ galaxies from the I. Labbé et al. (2017, in preparation) UltraVISTA selection in Appendix B—since they play a role in our volume density determination—such that this paper includes the full set of properties for the six most luminous $z \sim 8$ sources identified over UltraVISTA. Full details on the sample assembly and analysis of the complete sample are presented in I. Labbé et al. (2017, in preparation).

The HST/WFC3 observations were performed in the $z_{098}$, $J_{125}$ and $H_{160}$ bands (Figure 2) in 2017 March, for a total of 1 orbit per source, reaching depths of $\sim 25.8, 25.0, 24.5$ mag, respectively. One source (UVISTA-Y-1) also benefits from the archival data of the program SUSHI (PI: N. Suzuki; PID: 14808), which provides coverage in the $Y_{105}$ and $JH_{140}$ bands to 25.6 and 25.5 mag, respectively ($5\sigma, 1''/2$ aperture diameter; Figure 1). Leveraging the new HST images, we reprocessed the

existing ground- and space-based data, extracting accurate flux measurements with the mophongo software (Labbé et al. 2006, 2010a, 2010b, 2013, 2015). In our analysis we also now included the recently released ground-based *g*, *r*, *i*, *z*, and *Y* data from the ultradeep layer of the HyperSuprimeCam survey (Aihara et al. 2017a, 2017b).

Our analysis confirms the photometric redshift of three sources (UVISTA-Y-1, UVISTA-Y-5, and UVISTA-Y-6) to be $8.0 < z_{photo} < 8.7$ (Figure 3). Their measured luminosity $M_{UV} \sim -22.3$ makes them perhaps the brightest, most reliable galaxies at $z \sim 8-9$ identified to date. The uniquely deep optical, NIR, and Spitzer/IRAC data available for these sources is the reason for our high confidence in their nature (Figure 5). However, our analysis also demonstrates that the remaining two sources (UVISTA-J-1 and UVISTA-J-2) are very likely lower-redshift interlopers, with nominal redshifts of $z_{phot} \sim 2$ (Figure 4).

The three $z \sim 8$ candidate LBGs are characterized by average UV continuum slopes $\beta = -1.74 \pm 0.35$, consistent with lower-redshift ($z \sim 7$), similarly bright samples of LBGs of Bowler et al. (2017), suggesting a differential evolution of $\beta$ for the most luminous galaxies compared to $L^*$ or sub-$L^*$ galaxies at early cosmic epochs. Our $\beta$ are bluer than the extrapolations of measurements for lower-luminosity LBGs from CANDELS data, although the large uncertainties make them consistent at $<1\sigma$ (e.g., Bouwens et al. 2014; Figure 6), preventing us from deriving any further conclusion on differential evolution.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Distribution with redshift of the absolute magnitude of the bright LBGs detected so far at $z \gtrsim 7$. We include measurements of Oesch et al. (2014), Bouwens et al. (2015), Roberts-Borsani et al. (2016), Calvi et al. (2016), Bowler et al. (2017), and Ono et al. (2017) as indicated by the legend. Our candidate LBGs are marked as red points with error bars and lie at the high-luminosity end of all candidate $z \sim 8$ LBGs to date. Given the depth and wavelength coverage of the observations available for our sources, our bright sample arguably constitutes the brightest and most reliable sample of $z \sim 8-9$ galaxies to date.</figcaption></figure>

For our bright source UVISTA-Y-1 we measured a size of $R_{e,circ} = 0.9 \pm 0.3$ kpc, consistent with sizes of similarly luminous LBGs at $z \sim 7-10$, and suggesting very mild evolution over the first 1.5 Gyr of cosmic time (Figure 7).

Finally, using the $1/V_{max}$ formalism of Avni & Bahcall (1980), we computed the volume density $\Phi$ associated to the three candidate $z \sim 8$ LBGs. We find $\Phi = 8.49_{-4.60}^{+8.23} \times 10^{-7}$ Mpc$^{-3}$ mag$^{-1}$ at $M_{UV} = -22.21 \pm 0.25$. This constitutes the first measurement of the number density of $M_{UV} < -22$ mag LBGs at $z \sim 8$ based on actual detection of sources. We also estimated the volume density of star-forming galaxies at $z \sim 8$, including the full sample of six $M_{UV} \lesssim -22$ mag galaxies identified in I. Labbé et al. (2017, in preparation), three of which are presented in Appendix B. Together with the three candidate $z \sim 8$ LBGs confirmed by our main analysis, they constitute a complete sample of $z \sim 8$ LBGs with $M_{UV} \lesssim -22$ mag. The measured volume density associated to this complete sample is $\Phi = 17.3_{-6.9}^{+10.3} \times 10^{-7}$ Mpc$^{-3}$ mag$^{-1}$ at $M_{UV} = -22.21 \pm 0.25$ mag.

Unfortunately, given the large statistical uncertainties, we cannot use our current constraints on the bright end of the LF to discriminate between a Schechter or double-power-law form. Improvements on this front could come either from the detection of sources at even higher luminosities, where the discrepancies between the Schechter and the double-power-law form are larger, or from increased samples of galaxies within the current luminosity ranges, reducing the statistical uncertainties on the LF measurements.

The luminous galaxies that are the focus of this study play a path-finding role in the build-up of the galaxy population at

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Rest-frame UV continuum slope ($\beta$) vs. UV absolute magnitude for the three $z \gtrsim 8$ LBGs (red points with errorbars). Individual measurements for the sample of $z \sim 7$ LBGs of Bowler et al. (2017) are indicated by the purple points. The blue and yellow points mark the bi-weight median measurements of Bouwens et al. (2014) for the $z \sim 7$ and $z \sim 8$ samples, respectively. The corresponding best-fit relation is shown by the solid blue and yellow lines. The dashed lines indicate the extrapolations of these relations to the luminosities probed in this work. Our measurements of $\beta$ are consistent with the extrapolations of the $z \sim 7-8$ relations.</figcaption></figure>

early times. This work not only highlights the value of wide-area surveys for finding such rare and important objects, but it also further stresses the importance of the high-resolution imaging provided by HST in the study of the galaxy populations at early cosmic epochs. Such imaging enables the refinement of the photometric redshifts and helps to identify interlopers, resulting in cleaner samples. Our results, however, are still based on photometric redshifts from broadband imaging. A more robust picture inevitably requires spectroscopic confirmation. To this aim, we have started spectroscopic follow-up with Keck/MOSFIRE and VLT/X-Shooter. The sensitivities at observed optical/NIR wavelengths make it challenging with the current instrumentation. Bright high-redshift objects like those analyzed in this work offer two possibilities: (1) they constitute prime candidates for future spectroscopic follow-up with the *James Webb Space Telescope (JWST)*; and (2) their brightness together with refined photometric redshifts suggests them as valid targets for current ALMA observations, possibly resulting in spectroscopic confirmation even before the start of operation of the *JWST*.

The authors would like to thank the referee for their careful reading of the manuscript and for providing constructive comments that helped improving the quality of the paper. M.S. would like to thank Pieter van Dokkum, Jorryt Matthee, and Corentin Schreiber for useful discussion. K.I.C. 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 program ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. 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 Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton
<figure><img src="image_10.png" /></figure>

**Figure 7.** Top panel: evolution with redshift of the circularized effective radius of LBGs with $M_{UV} \sim -22$ mag. The size we measured for the bright $z \sim 8$ galaxy UVISTA-Y-1 is marked by the red filled circle. Also shown are measurements of the $R_{e,circ}$ for LBGs with $M_{UV} \sim -22$ from the literature, as indicated by the legend. To improve readability, the points of Huang et al. (2013) and of Curtis-Lake et al. (2016) have been arbitrarily shifted by $\Delta z = 0.1$ and $\Delta z = -0.1$, respectively. The dashed curve corresponds to the best fit of Shibuya et al. (2015) on the average measurements ($R_{e,circ} \propto (1+z)^{-0.84}$ kpc). Given the mild evolution in the characteristic luminosity of the UV LF of LBGs over $z \sim 4-10$ and the fact that we are considering the brightest end of the UV LF, the evolution presented here could roughly corresponds to the evolution of the descendants of objects like GN-z11, indicated by the filled triangle at $z \sim 11$. Bottom panel: evolution with redshift of the SFR surface density ($\Sigma_{SFR}$) for the same sample presented in the top panel, and computed following Ono et al. (2013). The data suggest a (slow) decrease of $\Sigma_{SFR}$ with cosmic time for the most luminous galaxies.

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 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 paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org. The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh,

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

**Figure 8.** Selection functions corresponding to the criteria adopted to identify candidate $z \sim 8$ LBGs, computed with a Monte Carlo simulation on real data. The top panel refers to the ultradeep stripes, and the lower panel to the UltraVISTA deep stripes. The original unevenly spaced measurements have been smoothed with a boxcar filter of 5 pixels, after casting them onto a regular grid for display purposes. Darker regions correspond to higher selection rates as indicated by the scale at the top of the figure. Our criteria allow us to select galaxies at $7.1 \lesssim z \lesssim 9.2$ and $6.9 \lesssim z \lesssim 9.3$, respectively for the ultradeep and deep regions.

Queens University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics, and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/IRFU, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Université de la Centre Nationale de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work is based in part on data products produced at Terapix available at the Canadian Astronomy Data Centre as part of the Canada-France-Hawaii Telescope Legacy Survey, a collaborative project of NRC and CNRS. This research has made
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Top panel: the red point with errorbars marks our estimate of volume density associated to the sample of candidate luminous $z \gtrsim 8$ galaxies considered in this work. The pink filled circle corresponds to the volume density after adding to the sample three luminous $z \sim 8$ LBGs from I. Labbé et al. (2017, in preparation), not targeted in the *HST* proposal. Recent UV LF determinations at $z \sim 8$ from Bradley et al. (2012), McLure et al. (2013), Schenker et al. (2013), Schmidt et al. (2014), Bouwens et al. (2015), and Finkelstein et al. (2015) are also reported, with plotting conventions detailed by the legend. The blue curve corresponds to the Schechter form of Bouwens et al. (2015), while the magenta curve represents the evolution to $z \sim 8$ of the $z \sim 7$ DPL of Bowler et al. (2017) according to Bouwens et al. (2016). Bottom panel: here we compare our volume density estimate (red point) to measurements of the UV LF at $z \sim 9$ (Oesch et al. 2013; Bouwens et al. 2016; Calvi et al. 2016; McLeod et al. 2016; Ishigaki et al. 2017). The blue curve represents the Schechter function from Bouwens et al. (2016). The magenta curve represents the bright end of the dual power law from Bowler et al. (2017) evolved to $z \sim 8$ following Bouwens et al. (2016) whose characteristic density has been adjusted to match that of the Schechter function at the characteristic luminosity.</figcaption></figure>

use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has benefited from the SpEx Prism Spectral Libraries, maintained by Adam Burgasser at http://pono.ucsd.edu/~adam/browndwarfs/spexprim.

## Appendix A
Inconsistent Flux Measurements and Interlopers

Our analysis contains three sources (UVISTA-Y-5, UVISTA-J-1, and UVISTA-J-2) which have flux measurements in the $H_{160}$ and UltraVISTA H bands that are inconsistent at $\gtrsim 2\sigma$. For UVISTA-J-1 and UVISTA-J-2, the lower flux measurement in the $H_{160}$ contributes to favoring a low-redshift solution. In this

section we discuss possible reasons for this systematic offset. Specifically we consider photometric zero-point offsets, time variability, nebular emission, and extended morphology.

When beginning the analysis for this paper, we had checked the zero point of the UltraVISTA DR3 H-band mosaic comparing the total flux in the $H_{160}$ and UltraVISTA H bands of $\sim 50$ bright unsaturated point sources, identified over the footprint of the CANDELS/COSMOS field. This region is covered by deep WFC3 data from the CANDELS program and by one of the ultradeep stripes of the UltraVISTA DR3 release, ensuring the highest S/N in both bands. The total flux was recovered from the curve-of-growth, measured out to $10''/0$ radii. In doing so, we masked any detected source, other than the point source itself, inside the radius adopted for the measurement. On this basis, we revised the UltraVISTA ZP estimates faintward by 0.1 mag. This zero-point offset is already included in the photometry we provide in Table 3.

In principle, the different releases of the UltraVISTA data set allow us to explore potential variations of source flux with time. We find that the photometry of UVISTA-J-1 in the UltraVISTA H band at the two epochs corresponding to DR1 and DR3 is consistent within uncertainties. Unfortunately, for UVISTA-Y-5 and UVISTA-J-2 this check is unfeasible, because the sources are undetected in the DR1 mosaic. While for UVISTA-Y-5 a non-detection in the DR1 data is still consistent at $\sim 2\sigma$ with the flux measured on the DR3 mosaic, this is not the case for UVISTA-J-2, whose flux should be $\sim 1.7 \times$ higher than UVISTA-Y-5, hence suggesting potential variability.

The coverage of the UltraVISTA H-band filter extends, at redder wavelengths, $\sim 10^3$ Å beyond that of the *HST*/WFC3 $H_{160}$. Red $H_{160}-H$ colors could then reflect the presence of nebular emission whose observed wavelength falls in the extra $\sim 10^3$ Å covered by the UltraVISTA H band. At $z \sim 8$ that wavelength range contains C III]λλ1907, 1909 Å. The resulting equivalent width would be in excess of $\sim 500$ Å, physically unlikely (see e.g., Stark et al. 2017 who found $EW_0 \sim 22$ Å for C III] at $z \sim 7$). Even assuming that C III] is dominating the flux in the UltraVISTA H band, similar EW for CIV would be required to explain the flux in the $K_S$ band, together with a very red SED red-ward of $K_S$ to match the fluxes in the 3.6 µm and 4.5 µm bands. One further possibility could be that the observed discrepancy was originated by [O III] + H β nebular lines if the galaxy had a redshift $z \sim 2.4-2.6$. In this case, $EW_0$([O III] + H β)$\gtrsim 1300$ Å. This value is larger than inferred from conversions of observed Hα EW at $z \sim 2$ (e.g., Erb et al. 2006); this solution becomes even more unlikely considering that the best-fit template required to match $z \sim 2$ is characterized by old age and dust attenuation. The old age would be inconsistent with strong nebular emission; furthermore, the dust attenuation would make the unreddened EW even larger. We conclude that nebular emission is not the likely origin of the observed discrepancy.

If the sources were characterized by extended, low surface-brightness wings, the short exposures in the WFC3 bands would not be enough to detect them, hence reducing their estimated luminosity. Unfortunately, a stacking of the WFC3 bands did not present any evidence for extended wings, possibly because of the low S/N characterizing our *HST* observations.

Visual inspection of the UltraVISTA mosaics showed that UVISTA-J-2 has a very bright ($\approx 12$ mag) neighboring star located $\sim 40$ arcsec north-east. Even though the procedure we adopted to extract the flux measurements takes care of
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Image stamps in inverted gray scale of the three additional bright candidate $z \gtrsim 8$ LBGs in the stacked optical, stacked Y, UltraVISTA J, H, and Ks, and Spitzer IRAC 3.6 and 4.5 µm bands. Each row corresponds to a source, as labeled on the left, where we omitted the prefix UVISTA from the object name for clarity. Each cutout is 5"0 × 5"0.</figcaption></figure>

<table><caption>**Table 6**
Total Flux Densities for the Three Additional Candidate $z \gtrsim 8$ LBGs over COSMOS/UltraVISTA Included in the Estimate of the LF</caption><thead><tr><th>Filter</th><th>UVISTA-Y-2<br>(nJy)</th><th>UVISTA-Y-3<br>(nJy)</th><th>UVISTA-Y-4<br>(nJy)</th></tr></thead><tbody><tr><td>CFHTLS u*</td><td>-6 ± 15</td><td>-6 ± 15</td><td>-10 ± 15</td></tr><tr><td>SSC B</td><td>-1 ± 10</td><td>-8 ± 8</td><td>5 ± 9</td></tr><tr><td>HSC g</td><td>2 ± 17</td><td>-11 ± 17</td><td>-8 ± 17</td></tr><tr><td>CFHTLS g</td><td>-1 ± 13</td><td>-12 ± 13</td><td>2 ± 13</td></tr><tr><td>SSC V</td><td>-14 ± 22</td><td>-14 ± 19</td><td>29 ± 20</td></tr><tr><td>HSC r</td><td>-4 ± 17</td><td>0 ± 16</td><td>7 ± 17</td></tr><tr><td>CFHTLS r</td><td>-12 ± 19</td><td>-4 ± 19</td><td>7 ± 20</td></tr><tr><td>SSC r+</td><td>-15 ± 20</td><td>-5 ± 17</td><td>10 ± 18</td></tr><tr><td>SSC i+</td><td>-139 ± 89</td><td>-77 ± 110</td><td>-103 ± 90</td></tr><tr><td>CFHTLS y</td><td>-18 ± 25</td><td>3 ± 24</td><td>-14 ± 25</td></tr><tr><td>CFHTLS i</td><td>-28 ± 26</td><td>-30 ± 25</td><td>3 ± 26</td></tr><tr><td>HSC i</td><td>-44 ± 24</td><td>-11 ± 24</td><td>-10 ± 24</td></tr><tr><td>CFHTLS z</td><td>-13 ± 55</td><td>-38 ± 54</td><td>7 ± 54</td></tr><tr><td>HSC z</td><td>-1 ± 35</td><td>-36 ± 35</td><td>-50 ± 35</td></tr><tr><td>SSC z+</td><td>-102 ± 75</td><td>-10 ± 73</td><td>43 ± 71</td></tr><tr><td>HSC y</td><td>31 ± 85</td><td>-89 ± 86</td><td>33 ± 84</td></tr><tr><td>UVISTA Y</td><td>29 ± 53</td><td>88 ± 53</td><td>173 ± 54</td></tr><tr><td>UVISTA J</td><td>410 ± 61</td><td>254 ± 61</td><td>432 ± 65</td></tr><tr><td>UVISTA H</td><td>432 ± 78</td><td>357 ± 77</td><td>392 ± 86</td></tr><tr><td>UVISTA K<sub>S</sub></td><td>275 ± 86</td><td>263 ± 84</td><td>266 ± 110</td></tr><tr><td>IRAC 3.6 µm</td><td>492 ± 50</td><td>589 ± 45</td><td>620 ± 68</td></tr><tr><td>IRAC 4.5 µm</td><td>799 ± 57</td><td>729 ± 49</td><td>682 ± 108</td></tr><tr><td>IRAC 5.8 µm</td><td>688 ± 1702</td><td>4546 ± 2069</td><td>-1686 ± 1819</td></tr><tr><td>IRAC 8.0 µm</td><td>1384 ± 2105</td><td>-2896 ± 2461</td><td>-795 ± 2123</td></tr></tbody></table>

**Note.** These measurements are reprocessed fluxes using HSC z band as prior for mophongo.

estimating the background, we cannot exclude a residual contamination from the wings of the bright source. It is remarkable, though, that if we exclude the UltraVISTA data from the SED of UVISTA-J-2, the *HST* photometry is consistent with the SED of an LBG at $z \sim 9.5$. Similarly, if we exclude the UltraVISTA data from UVISTA-J-1, we obtain an SED consistent with an LBG at $z \sim 9.5$, although less

robust than UVISTA-J-2 given the detection in the HSC z band. In conclusion this suggests that for these two objects a high-redshift solution is not completely ruled out.

As we show in Figure 4, the WFC3 observations are responsible for, or contribute to, the solution at $z \sim 2$. The above considerations stress the importance of performing high-S/N *HST* follow-up of the candidate bright LBGs detected from ground-based surveys, in order to reduce the systematic uncertainties in the photometry and produce more stable photometric redshift measurements. The current samples of $z \gtrsim 4$ LBGs at lower luminosities are generally built from deep *HST* imaging. The higher S/N of the *HST* observations in these fields greatly reduces the chances of uncertain redshifts identification of their redshifts. Nonetheless, issues in the assessment of the nature of candidate LBGs arise at the faint end of the UV LF (see e.g., Bouwens et al. 2017).

## Appendix B
Sources Used to Estimate the LF Not Included in Our *HST* Follow-up Program

Here we present the three additional candidate bright $z \sim 8$ LBGs, from the sample of I. Labbé et al. (2017, in preparation), that we included in our LF estimates (Section 5.4). Their selection followed the same methods as described in Section 2. However, due to the lack of *HST* imaging, we reprocessed their photometry with *mophongo* adopting the HSC z band as positional and morphological prior: its red effective wavelength together with its depth allowed us to detect almost every source (i.e., neighboring, potentially contaminating objects) on the UltraVISTA and IRAC mosaics, while its narrow PSF (the narrowest among the ground-based and IRAC data sets) ensures we can consistently use it as prior with *mophongo* to recover the flux in all the bands.

The fact that our targets are *Y* or *J* dropouts, undetected by construction in the z band, does not constitute a major problem for our photometry. Indeed, *mophongo* can perform the aperture photometry blindly, i.e., without the need to detect the source of interest.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Left panels: SEDs for the three additional candidate LBGs at $z \gtrsim 8$, included in the estimate of the LF. The colored squares with black errorbars mark the photometric measurements, while arrows represent $2\sigma$ upper limits. Open squares and arrows mark the IRAC 3.6, 4.5, and 5.8 µm bands, not used for the measurement of the fiducial photometric redshift. Photometry in the HSC Survey bands is represented in yellow. The fiducial best-fit SED template from EAZY is indicated by the thick blue curve; the thin dark blue curve represents the best-fit SED when all bands are used for the photometric redshift measurement. Also shown are the best-fitting brown dwarf template (light brown curve) and the solution when the redshift is forced to be $z < 6$ (gray curve). Right panels: redshift likelihood distributions ($p(z)$) for the three LBGs for the fiducial solution (blue) and for the solution obtained considering the full set of flux measurements. The label in the top-left corner indicates the estimated photometric redshifts. The $p(z)$ are peaked, with no integrated probability for a secondary solution at lower redshifts.</figcaption></figure>

<table><caption>**Table 7**
Main Properties of the Three Additional Bright Candidate $z \sim 8$ LBGs Included in the LF Measurement</caption><thead><tr><th>ID</th><th>R.A.<br>(J2000)</th><th>Decl.<br>(J2000)</th><th>H<br>(mag)</th><th>z<sub>photo</sub><sup>a</sup></th><th>p(z > 6)<sup>b</sup></th><th>M<sub>UV</sub><sup>c</sup><br>(mag)</th><th>β<sup>d</sup></th><th>EW<sub>0</sub>(Hβ + [O III])<sup>e</sup><br>(Å)</th></tr></thead><tbody><tr><td>UVISTA-Y-2</td><td>10:02:12.558</td><td>2:30:45.71</td><td>24.8</td><td>8.30<sup>+0.47</sup><sub>-0.50</sub></td><td>1.00</td><td>-22.37 ± 0.20</td><td>-2.45 ± 0.57</td><td>705<sup>+352</sup><sub>-287</sub></td></tr><tr><td>UVISTA-Y-3</td><td>10:00:32.322</td><td>1:44:31.26</td><td>25.0</td><td>8.24<sup>+0.81</sup><sub>-0.86</sub></td><td>1.00</td><td>-22.02 ± 0.23</td><td>-1.93 ± 0.74</td><td>269<sup>+219</sup><sub>-188</sub></td></tr><tr><td>UVISTA-Y-4</td><td>10:00:58.485</td><td>1:49:55.96</td><td>24.9</td><td>7.42<sup>+0.20</sup><sub>-0.23</sub></td><td>1.00</td><td>-22.18 ± 0.24</td><td>-2.69 ± 0.68</td><td>124<sup>+442</sup><sub>-331</sub></td></tr></tbody></table>

**Notes.**

In Figure 10 we present image cutouts in the optical, NIR, and IRAC 3.6 and 4.5 µm bands. Table 6 presents the flux density measurements of these three objects, while Table 7 lists

their main observational properties obtained following the same analysis adopted for our main sample. Their observed and best-fit SEDs and $p(z)$ are shown in Figure 11. The $\chi^2$ for the

<footnote><sup>a</sup> Best photometric redshift estimate from EAZY, excluding the IRAC bands from the fit, and corresponding 68% confidence interval.</footnote>

<footnote><sup>b</sup> Probability, computed by EAZY, that the solution is at $z > 6$.</footnote>

<footnote><sup>c</sup> Absolute magnitudes at rest-frame 1600 Å from EAZY.</footnote>

<footnote><sup>d</sup> Rest-frame UV continuum slopes from the UltraVISTA *J*, *H*, and *K*<sub>S</sub> bands.</footnote>

<footnote><sup>e</sup> Rest-frame equivalent width of Hβ + [O III] obtained from the [3.6] to [4.5] color assuming an SED flat in $f_ν$ (i.e., $β = -2$).</footnote>
$z \sim 8$ solutions are $\chi^2 = 17.6, 28.7, 16.6$, for UVISTA-Y-2, UVISTA-Y-3, and UVISTA-Y-4, respectively. The sources are characterized by blue UV continuum slopes ($\beta \sim -2$), excluding red/dusty interloper solutions. Indeed, forcing $z < 6$ generates solutions with $\chi^2 = 76.7, 56.7, 50.1$. Our measured photometry allows us to exclude a solution as brown dwarves as well ($\chi^2 = 70.6, 76.5, 57.0)$.

Their relative brightness ($H \sim 25$ mag) then translates into high UV luminosities, with absolute magnitudes $M_{UV} \lesssim -22$ mag.

## ORCID iDs

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

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

Gabriel B. Brammer https://orcid.org/0000-0003-2680-005X

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

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

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

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

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

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

Benne W. Holwerda https://orcid.org/0000-0002-4884-6756

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

Renske Smit https://orcid.org/0000-0001-8034-7802

## References

Aihara, H., Arimoto, N., Armstrong, R., et al. 2017a, PASJ, in press (arXiv:1704.05858)

Aihara, H., Armstrong, R., Bickerton, S., et al. 2017b, PASJ, in press (arXiv:1702.08449)

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

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

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

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

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

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

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., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017, ApJ, 843, 129

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

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

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

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

Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13

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

Burgasser, A. J. 2014, in Ast. Soc. India Conf. Ser, 11, 7

Burrows, A., Sudarsky, D., & Hubeny, I. 2006, ApJ, 640, 1063

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

Cappelluti, N., Brusa, M., Hasinger, G., et al. 2009, A&A, 497, 635

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

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

Curtis-Lake, E., McLure, R. J., Dunlop, J. S., et al. 2016, MNRAS, 457, 440

Duncan, K., & Conselice, C. J. 2015, MNRAS, 451, 2030

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Dunlop, J. S., Rogers, A. B., McLure, R. J., et al. 2015, MNRAS, 432, 3520

Erb, D. K., Steidel, C. C., Shapley, A. E., et al. 2006, ApJ, 647, 128

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

Fernández-Soto, A., Lanzetta, K. M., & Yahil, A. 1999, ApJ, 513, 34

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

Gehrels, N. 1986, ApJ, 303, 336

Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35

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

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

Holwerda, B. W., Trenti, M., Clarkson, W., et al. 2014, ApJ, 788, 77

Huang, K.-H., Ferguson, H. C., Ravindranath, S., & Su, J. 2013, ApJ, 765, 68

Huertas-Company, M., Pérez-González, P. G., Mei, S., et al. 2015, ApJ, 809, 95

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

Ishigaki, M., Kawamata, R., Ouchi, M., Oguri, M., & Shimasaku, K. 2017, arXiv:1702.04867

Jiang, L., Egami, E., Mechtley, M., et al. 2013, ApJ, 772, 99

Kennicutt, R. C., Jr. 1998, ARA&A, 36, 189

Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36

Koekemoer, A. M., Fruchter, A. S., Hook, R. N., & Hack, W. 2003, in Proc. HST Calibration Workshop. Hubble after the Installation of the ACS and the NICMOS Cooling System, ed. S. Arribas, A. Koekemoer, & B. Whitmore (Baltimore, MD: STScI), 337

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

Laidler, V. G., Papovich, C., Grogin, N. A., et al. 2007, PASP, 119, 1325

Law, D. R., Steidel, C. C., Erb, D. K., et al. 2007, ApJ, 656, 1

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

Marchesi, S., Civano, F., Elvis, M., et al. 2016, ApJ, 817, 34

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

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

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

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

Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15

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

Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185

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

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

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

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

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

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

Ono, Y., Ouchi, M., Harikane, Y., et al. 2017, PASJ, in press (arXiv:1704.06004)

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097

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

Ribeiro, B., Le Fèvre, O., Tasca, L. A. M., et al. 2016, A&A, 593, A22

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

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

Ryan, R. E., Thorman, P. A., Yan, H., et al. 2011, ApJ, 739, 83

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

Schechter, P. 1976, ApJ, 203, 297

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

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

Schmidt, M. 1968, ApJ, 151, 393

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

Sërsic, J. L. 1968, Atlas de Galaxias Australes (Cordoba: Observatorio Astronomico)

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

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

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

Smolcic, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

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., Schenker, M. A., Ellis, R., et al. 2013, ApJ, 763, 129

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2017, ApJ, 843, 86

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

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

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

van Dokkum, P. G., Brammer, G., Fumagalli, M., et al. 2011, ApJL, 743, L15

Wilkins, S. M., Bouwens, R. J., Oesch, P. A., et al. 2016, MNRAS, 455, 659

Wilkins, S. M., Bunker, A. J., Stanway, E., Lorenzoni, S., & Caruana, J. 2011, MNRAS, 417, 717

Williams, R. E., Blacker, B., Dickinson, M., et al. 1996, AJ, 112, 1335

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