HST imaging of the brightest z~8-9 galaxies from UltraVISTA: the extreme bright end of the UV luminosity function

TL;DR

This study discovers and confirms three bright $z extasciitilde8$-9 galaxies using HST imaging, measuring their properties and estimating the volume density at the bright end of the early universe's UV luminosity function.

Contribution

It provides new measurements of bright $z extasciitilde8$-9 galaxies and estimates their volume density, advancing understanding of the UV luminosity function at high redshift.

Findings

Confirmed 3 galaxies at $z extasciitilde8$-9 with HST data.

Measured mean UV slopes $eta extasciitilde -1.91$ for these galaxies.

Estimated the volume density of extremely bright galaxies at $z extasciitilde8$.

Abstract

We report on the discovery of three especially bright candidate $z_{phot} \gtrsim 8$ galaxies. Five sources were targeted for follow-up with HST/WFC3, selected from a larger sample of 16 bright ($24.8 \lesssim H\lesssim25.5$~mag) candidate $z\gtrsim 8$ LBGs identified over the 1.6 degrees$^2$ of the COSMOS/UltraVISTA field. These were identified as Y and J dropouts by leveraging the deep (Y-to-$K_{S} \sim 25.3-24.8$~mag, $5\sigma$) NIR data from the UltraVISTA DR3 release, deep ground based optical imaging from the CFHTLS and Subaru 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 HyperSuprime Cam g, r, i, z and Y band data, we confirm that 3/5 galaxies have robust $z_{phot}\sim8.0-8.7$, consistent with the initial selection. The remaining 2/5 galaxies have a nominal $z_{phot}\sim2$. However, if we use the HST data alone, these objects have increased probability of being at $z\sim9$. Furthermore, we measure mean UV continuum slopes $\beta=-1.91\pm0.26$ for the three $z\sim8-9$ galaxies, marginally bluer than similarly luminous $z\sim4-6$ in CANDELS but consistent with previous measurements of similarly luminous galaxies at $z\sim7$. The circularized effective radius for our brightest source is $0.9\pm0.2$ kpc, similar to previous measurements for a bright $z\sim11$ galaxy and bright $z\sim7$ galaxies. Finally, enlarging our sample to include the six brightest $z\sim8$ LBGs identified over UltraVISTA (i.e., including three other sources from Labbe et al. 2017, in prep.) we estimate for the first time the volume density of galaxies at the extreme bright ($M_{UV}\sim-22$~mag) end of the $z\sim8$ UV LF. Despite this exceptional result, the still large statistical uncertainties do not allow us to discriminate between a Schechter and a double power-law form.

Tables

Table 1. Table 1: Photometric depths of the adopted ground-based and Spitzer /IRAC data sets, and corresponding average aperture corrections.

Filter	Aperture	Depth
name	correctionaaAverage multiplicative factors applied to estimate total fluxes.	$5\sigma$bbAverage depth over the full field corresponding to $5\sigma$ flux dispersions in empty apertures of $1\stackrel{\prime \prime }{.}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\mu$m and $4.5\mu$m bands correspond to the regions with SMUVS+SCOSMOS+SPLASH coverage (approximately overlapping with the ultradeep stripes) and SPLASH+SCOSMOS only ($\approx$ deep stripes).
CFHTLS $u^{\ast }$	$2.3$	$26.7$
SSC $B$	$1.7$	$27.4$
HSC $g$ccThe HyperSuprimeCam data was 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.	$2.1$	$26.7$
CFHTLS $g$	$2.2$	$26.8$
SSC $V$	$2.2$	$26.4$
HSC $r$ccThe HyperSuprimeCam data was 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.	$1.6$	$26.8$
CFHTLS $r$	$2.1$	$26.4$
SSC $r^{+}$	$2.0$	$26.6$
SSC $i^{+}$	$1.9$	$26.2$
CFHTLS $y$	$2.0$	$26.1$
CFHTLS $i$	$2.0$	$26.0$
HSC $i$ccThe HyperSuprimeCam data was 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.	$1.7$	$26.3$
CFHTLS $z$	$2.1$	$25.2$
HSC $z$ccThe HyperSuprimeCam data was 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.	$1.6$	$25.9$
SSC $z^{+}$	$2.3$	$25.0$
HSC $y$ccThe HyperSuprimeCam data was 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.	$2.1$	$24.9$
UVISTA $Y$	$2.6$	$25.4/24.5$
UVISTA $J$	$2.4$	$25.4/24.4$
UVISTA $H$	$2.2$	$25.1/24.1$
UVISTA $K_S$	$2.2$	$24.8/23.7$
IRAC $3.6\mu$m	$5.3$	$24.9/24.5$
IRAC $4.5\mu$m	$5.4$	$24.7/24.3$
IRAC $5.8\mu$m	$8.4$	$20.8$
IRAC $8.0\mu$m	$10.1$	$20.6$

Table 2. Table 2: HST observations we obtained over the bright z ∼ 8 − 9 similar-to 𝑧 8 9 z\sim 8-9 candidates from Labbé et al. ( 2017, in preparation ) .

ID	R.A.	Dec.	$H$	Exposure times					Photometric depths
				$z_{098}$	$Y_{105}$aaFortuitously, observations in the $Y_{105}$ and $J{H}_{140}$ bands are available over over one $z\sim 8-9$ candidate in our program as part of the separate HST program SUSHI (PI: Nao Suzuki - PID: 14808).	$J_{125}$	$J{H}_{140}$aaFortuitously, observations in the $Y_{105}$ and $J{H}_{140}$ bands are available over over one $z\sim 8-9$ candidate in our program as part of the separate HST program SUSHI (PI: Nao Suzuki - PID: 14808).	$H_{160}$	$z_{098}$	$Y_{105}$aaFortuitously, observations in the $Y_{105}$ and $J{H}_{140}$ bands are available over over one $z\sim 8-9$ candidate in our program as part of the separate HST program SUSHI (PI: Nao Suzuki - PID: 14808).	$J_{125}$	$J{H}_{140}$aaFortuitously, observations in the $Y_{105}$ and $J{H}_{140}$ bands are available over over one $z\sim 8-9$ candidate in our program as part of the separate HST program SUSHI (PI: Nao Suzuki - PID: 14808).	$H_{160}$
	[J2000]	[J2000]	[mag]	[sec]	[sec]	[sec]	[sec]	[sec]	[mag]	[mag]	[mag]	[mag]	[mag]
UVISTA-Y-1	09:57:47.90	2:20:43.7	$24.8$	$1512$	1022	$462$	$1197$	$412$	$25.8$	$25.6$	$25.0$	$25.5$	$24.6$
UVISTA-Y-5	10:00:31.89	1:57:50.2	$24.9$	$1512$	$\cdots$	$462$	$\cdots$	$412$	$25.8$	$\cdots$	$25.0$	$\cdots$	$24.6$
UVISTA-Y-6	10:00:12.51	2:03:00.5	$25.3$	$1512$	$\cdots$	$462$	$\cdots$	$412$	$25.8$	$\cdots$	$25.0$	$\cdots$	$24.6$
UVISTA-J-1	10:02:25.48	2:29:13.6	$24.6$	$1512$	$\cdots$	$462$	$\cdots$	$412$	$25.8$	$\cdots$	$25.0$	$\cdots$	$24.6$
UVISTA-J-2	09:59:07.19	1:56:54.0	$24.4$	$1512$	$\cdots$	$462$	$\cdots$	$412$	$25.8$	$\cdots$	$25.0$	$\cdots$	$24.6$

Table 3. Table 3: Total flux densities for the five candidate z ≳ 8 greater-than-or-equivalent-to 𝑧 8 z\gtrsim 8 LBGs over COSMOS/UltraVISTA targeted by our small HST program.

Filter	UVISTA-Y-1		UVISTA-Y-5		UVISTA-Y-6		UVISTA-J-1		UVISTA-J-2
	[nJy]		[nJy]		[nJy]		[nJy]		[nJy]
CFHTLS $u^{\ast }$	$\cdots$		$-13$ $\pm$	$17$	$8$ $\pm$	$14$	$4$ $\pm$	$15$	$13$ $\pm$	$17$
SSC $B$	$-2$ $\pm$	$7$	$4$ $\pm$	$9$	$-11$ $\pm$	$9$	$-4$ $\pm$	$8$	$5$ $\pm$	$11$
HSC $g$	$8$ $\pm$	$15$	$-10$ $\pm$	$20$	$4$ $\pm$	$16$	$1$ $\pm$	$15$	$6$ $\pm$	$20$
CFHTLS $g$	$\cdots$		$3$ $\pm$	$15$	$1$ $\pm$	$13$	$-8$ $\pm$	$10$	$-4$ $\pm$	$17$
SSC $V$	$-12$ $\pm$	$17$	$-1$ $\pm$	$24$	$0$ $\pm$	$21$	$-7$ $\pm$	$17$	$-11$ $\pm$	$27$
HSC $r$	$-7$ $\pm$	$14$	$-3$ $\pm$	$18$	$11$ $\pm$	$15$	$4$ $\pm$	$15$	$10$ $\pm$	$17$
CFHTLS $r$	$\cdots$		$-15$ $\pm$	$23$	$8$ $\pm$	$21$	$-22$ $\pm$	$17$	$7$ $\pm$	$22$
SSC $r^{\prime }$	$7$ $\pm$	$16$	$-29$ $\pm$	$23$	$20$ $\pm$	$19$	$15$ $\pm$	$16$	$-7$ $\pm$	$18$
CFHTLS $y$	$\cdots$		$-11$ $\pm$	$28$	$22$ $\pm$	$26$	$16$ $\pm$	$23$	$-26$ $\pm$	$28$
CFHTLS $i$	$\cdots$		$-11$ $\pm$	$29$	$5$ $\pm$	$29$	$-28$ $\pm$	$24$	$-29$ $\pm$	$29$
HSC $i$	$21$ $\pm$	$21$	$-20$ $\pm$	$27$	$1$ $\pm$	$23$	$19$ $\pm$	$21$	$23$ $\pm$	$27$
SSC $i^{+}$	$-9$ $\pm$	$22$	$-36$ $\pm$	$26$	$2$ $\pm$	$21$	$-18$ $\pm$	$24$	$-24$ $\pm$	$27$
CFHTLS $z$	$\cdots$		$6$ $\pm$	$63$	$-13$ $\pm$	$59$	$-65$ $\pm$	$52$	$32$ $\pm$	$62$
HSC $z$	$9$ $\pm$	$31$	$-27$ $\pm$	$39$	$17$ $\pm$	$33$	$52$ $\pm$	$31$	$8$ $\pm$	$37$
SSC $z^{\prime }$	$51$ $\pm$	$64$	$-51$ $\pm$	$93$	$69$ $\pm$	$85$	$70$ $\pm$	$64$	$-39$ $\pm$	$83$
HSC $Y$	$-31$ $\pm$	$76$	$-91$ $\pm$	$98$	$89$ $\pm$	$80$	$80$ $\pm$	$76$	$-29$ $\pm$	$98$
$z_{098}$	$46$ $\pm$	$34$	$-7$ $\pm$	$34$	$9$ $\pm$	$34$	$13$ $\pm$	$34$	$22$ $\pm$	$34$
UVISTA $Y$	$18$ $\pm$	$48$	$-42$ $\pm$	$68$	$16$ $\pm$	$51$	$67$ $\pm$	$55$	$37$ $\pm$	$67$
$Y_{105}$	$92$ $\pm$	$41$	$\cdots$		$\cdots$		$\cdots$		$\cdots$
$J_{125}$	$279$ $\pm$	$70$	$195$ $\pm$	$70$	$172$ $\pm$	$70$	$220$ $\pm$	$70$	$304$ $\pm$	$70$
UVISTA $J$	$324$ $\pm$	$50$	$235$ $\pm$	$66$	$211$ $\pm$	$53$	$125$ $\pm$	$50$	$195$ $\pm$	$66$
$J{H}_{140}$	$303$ $\pm$	$44$	$\cdots$		$\cdots$		$\cdots$		$\cdots$
$H_{160}$	$511$ $\pm$	$107$	$152$ $\pm$	$107$	$272$ $\pm$	$107$	$188$ $\pm$	$107$	$257$ $\pm$	$107$
UVISTA $H$	$455$ $\pm$	$61$	$393$ $\pm$	$86$	$280$ $\pm$	$66$	$510$ $\pm$	$70$	$657$ $\pm$	$86$
UVISTA $K_S$	$480$ $\pm$	$77$	$321$ $\pm$	$102$	$271$ $\pm$	$82$	$602$ $\pm$	$95$	$822$ $\pm$	$102$
IRAC $3.6\mu$m	$623$ $\pm$	$85$	$289$ $\pm$	$74$	$434$ $\pm$	$106$	$1162$ $\pm$	$90$	$913$ $\pm$	$110$
IRAC $4.5\mu$m	$931$ $\pm$	$109$	$589$ $\pm$	$86$	$598$ $\pm$	$130$	$1277$ $\pm$	$102$	$1204$ $\pm$	$130$
IRAC $5.8\mu$m	$-2893$ $\pm$	$2568$	$-1978$ $\pm$	$4831$	$-643$ $\pm$	$3000$	$-2209$ $\pm$	$2655$	$8075$ $\pm$	$4602$
IRAC $8.0\mu$m	$1423$ $\pm$	$3021$	$499$ $\pm$	$6310$	$-3325$ $\pm$	$3803$	$1992$ $\pm$	$3771$	$6734$ $\pm$	$5613$

Table 4. Table 4: Physical parameters for the three galaxies with photometric redshift z ≳ 8 greater-than-or-equivalent-to 𝑧 8 z\gtrsim 8 $\mathsection$ $\mathsection$ footnotemark:

ID	$z_{\mathrm{photo}}$aaBest photometric redshift estimate from EAzY, excluding the IRAC bands from the fit, and corresponding 68% confidence interval.	$p\left(z>6\right)$   bbProbability, computed by EAzY, that the solution is at $z>6$.	$M_{\mathrm{UV}}$ccAbsolute magnitudes at rest-frame 1600Å from EAzY.	$\beta$ddRest-frame UV continuum slopes from the $H_{160}$ and UltraVISTA $H$ and $K_S$ bands.	EW${}_0{}^{}(H\beta +\left[\mathrm{OIII}\right])$eeRest-frame equivalent width of H$\beta +\left[\mathrm{OIII}\right]$ obtained from the $K_S-\left[4.5\right]$ color assuming an SED flat in $f_{\nu }$ (i.e. $\beta =-2$). If $\beta =-1.8$, the EW and associated uncertainties would be a factor $\sim 1.2$ smaller.
			[mag]		[Å]
UVISTA-Y-1	${8.02}_{-0.49}^{+0.41}$	$0.88$	$-22.46\pm 0.15$	$-1.98\pm 0.67$	${1041}_{-515}^{+713}$
UVISTA-Y-5$\dagger$$\dagger$footnotemark:	${8.39}_{-0.60}^{+0.60}$	$0.99$	$-22.37\pm 0.24$	$-1.62\pm 0.42$	${887}_{-686}^{+1323}$
UVISTA-Y-6	${8.35}_{-0.81}^{+0.65}$	$0.83$	$-21.97\pm 0.26$	$-2.12\pm 1.44$	${1291}_{-940}^{+1749}$

Table 5. Table 5: Morphological parameters for UVISTA-Y-1 measured on the J ​ H 140 𝐽 subscript 𝐻 140 JH_{140} band with galfit and SExtractor

Algorithm	R.A.	Dec.	$R_{e,\mathrm{circ}}$aaCircularized effective radius	$q$bbMinor-to-major axis ratio	$n$ccSérsic index. This was kept fixed when running galfit	${\Sigma }_{\mathrm{SFR}}$ddStar-formation rate surface density, computed following Ono et al. (2013)
	[J2000]	[J2000]	[kpc]			$[M_{\odot }$ yr-1 kpc-2]
galfit	09:57:47.910	+02:20:43.50	$0.9\pm 0.3$	$0.9\pm 0.2$	1.5	${11}_{-5}^{+14}$
SExtractor	09:57:47.910	+02:20:43.50	$0.7\pm 0.1$	$0.6\pm 0.1$	$\cdots$	${17}_{-5}^{+10}$

Table 6. Table 6: Total flux densities for the three additional candidate z ≳ 8 greater-than-or-equivalent-to 𝑧 8 z\gtrsim 8 LBGs over COSMOS/UltraVISTA included in the estimate of the LF.

Filter	UVISTA-Y-2		UVISTA-Y-3		UVISTA-Y-4
	[nJy]		[nJy]		[nJy]
CFHTLS $u^{\ast }$	$-6$ $\pm$	$15$	$-6$ $\pm$	$15$	$-10$ $\pm$	$15$
SSC $B$	$-1$ $\pm$	$10$	$-8$ $\pm$	$8$	$5$ $\pm$	$9$
HSC $g$	$2$ $\pm$	$17$	$-11$ $\pm$	$17$	$-8$ $\pm$	$17$
CFHTLS $g$	$-1$ $\pm$	$13$	$-12$ $\pm$	$13$	$2$ $\pm$	$13$
SSC $V$	$-14$ $\pm$	$22$	$-14$ $\pm$	$19$	$29$ $\pm$	$20$
HSC $r$	$-4$ $\pm$	$17$	$0$ $\pm$	$16$	$7$ $\pm$	$17$
CFHTLS $r$	$-12$ $\pm$	$19$	$-4$ $\pm$	$19$	$7$ $\pm$	$20$
SSC $r^{+}$	$-15$ $\pm$	$20$	$-5$ $\pm$	$17$	$10$ $\pm$	$18$
SSC $i^{+}$	$-139$ $\pm$	$89$	$-77$ $\pm$	$110$	$-103$ $\pm$	$90$
CFHTLS $y$	$-18$ $\pm$	$25$	$3$ $\pm$	$24$	$-14$ $\pm$	$25$
CFHTLS $i$	$-28$ $\pm$	$26$	$-30$ $\pm$	$25$	$3$ $\pm$	$26$
HSC $i$	$-44$ $\pm$	$24$	$-11$ $\pm$	$24$	$-10$ $\pm$	$24$
CFHTLS $z$	$-13$ $\pm$	$55$	$-38$ $\pm$	$54$	$7$ $\pm$	$54$
HSC $z$	$-1$ $\pm$	$35$	$-36$ $\pm$	$35$	$-50$ $\pm$	$35$
SSC $z^{+}$	$-102$ $\pm$	$75$	$-10$ $\pm$	$73$	$43$ $\pm$	$71$
HSC $y$	$31$ $\pm$	$85$	$-89$ $\pm$	$86$	$33$ $\pm$	$84$
UVISTA $Y$	$29$ $\pm$	$53$	$88$ $\pm$	$53$	$173$ $\pm$	$54$
UVISTA $J$	$410$ $\pm$	$61$	$254$ $\pm$	$61$	$432$ $\pm$	$65$
UVISTA $H$	$432$ $\pm$	$78$	$357$ $\pm$	$77$	$392$ $\pm$	$86$
UVISTA $K_S$	$275$ $\pm$	$86$	$263$ $\pm$	$84$	$266$ $\pm$	$110$
IRAC $3.6\mu$m	$492$ $\pm$	$50$	$589$ $\pm$	$45$	$620$ $\pm$	$68$
IRAC $4.5\mu$m	$799$ $\pm$	$57$	$729$ $\pm$	$49$	$682$ $\pm$	$108$
IRAC $5.8\mu$m	$688$ $\pm$	$1702$	$4546$ $\pm$	$2069$	$-1686$ $\pm$	$1819$
IRAC $8.0\mu$m	$1384$ $\pm$	$2105$	$-2896$ $\pm$	$2461$	$-795$ $\pm$	$2123$

Table 7. Table 7: Main properties of the three additional bright candidate z ∼ 8 similar-to 𝑧 8 z\sim 8 LBG included in the LF measurement.

ID	R.A.	Dec.	H	$z_{\mathrm{photo}}$aaBest photometric redshift estimate from EAzY, excluding the IRAC bands from the fit, and corresponding 68% confidence interval.	$p\left(z>6\right)$   bbProbability, computed by EAzY, that the solution is at $z>6$.	$M_{\mathrm{UV}}$ccAbsolute magnitudes at rest-frame 1600Å from EAzY.	$\beta$ddRest-frame UV continuum slopes from the UltraVISTA $J$, $H$ and $K_S$ bands.	$E{W}_0\left(H\beta +\left[\mathrm{OIII}\right]\right)$eeRest-frame equivalent width of H$\beta +\left[\mathrm{OIII}\right]$ obtained from the $\left[3.6\right]-\left[4.5\right]$ color assuming an SED flat in $f_{\nu }$ (i.e. $\beta =-2$).
	[J2000]	[J2000]	[mag]			[mag]		[Å]
UVISTA-Y-2	10:02:12.558	2:30:45.71	$24.8$	${8.30}_{-0.50}^{+0.47}$	$1.00$	$-22.37\pm 0.20$	$-2.45\pm 0.57$	${705}_{-287}^{+352}$
UVISTA-Y-3	10:00:32.322	1:44:31.26	$25.0$	${8.24}_{-0.86}^{+0.81}$	$1.00$	$-22.02\pm 0.23$	$-1.93\pm 0.74$	${269}_{-188}^{+219}$
UVISTA-Y-4	10:00:58.485	1:49:55.96	$24.9$	${7.42}_{-0.23}^{+0.20}$	$1.00$	$-22.18\pm 0.24$	$-2.69\pm 0.68$	${124}_{-331}^{+412}$

Equations

$$Y-(J+H)/2>0.75\vee (J-H)>0.8$$

$$(H-K_{\mathrm{S}}<0.7)\wedge (K_{\mathrm{S}}-[3.6]<1.75\vee H-[3.6]<1.75)$$

Full text

galaxies: formation, galaxies: evolution, 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 $\sim 1500$ candidate galaxies at $z>6$, $\sim 200$ of which at $z\sim 8-10$ (e.g., Oesch et al. 2013; Bouwens et al. 2015; Finkelstein et al. 2015; Bouwens et al. 2016b). These samples are characterised by $M_{\mathrm{UV}}\gtrsim-22$ (i.e., apparent magnitudes fainter than $\gtrsim 25.5$ at $z\sim 8$), $\sim 1$ mag more luminous than the current determinations of the characteristic magnitude of the rest-frame ultraviolet (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; McLure et al. 2013; Schenker et al. 2013; Duncan et al. 2014; Bouwens et al. 2015; Finkelstein et al. 2015), galaxies fainter than the characteristic luminosity dominate the estimates of the star-formation rate 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_{\mathrm{UV}}=-22.1$ mag) galaxy at the record redshift of $z_{\mathrm{grism}}=11.09_{-0.12}^{+0.08}$, 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_{\mathrm{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 deg2, 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_{\mathrm{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 to 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, b). Analysis of this three-layered dataset resulted in a sample of $\sim 600$ $z\sim 6-7$ LBG galaxies ($\sim 70$ galaxies at $z\sim 7$) with $M_{\mathrm{UV}}\lesssim-25$ mag over $\sim 100$ deg2. After carefully removing 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 CFHTLS (Erben et al., 2009; Hildebrandt et al., 2009) and Subaru Suprime-Cam (Taniguchi et al. 2007), and with deep IRAC mosaics 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$ (Labbé et al., 2017, in preparation).

The primary question is if 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_{\mathrm{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 under-estimate 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 suppressing low signal-to-noise (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 dataset, 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 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) in order 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_{\mathrm{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 Sect. 2 we briefly describe the datasets and the selection criteria adopted for the assembly of the sample; in Sect. 3 we describe the HST data and how the photometry was performed; the results are presented in Sect. 4; in Sect. 5 we discuss the results and in Sect. 6 we conclude.

Throughout this work, we use the following short-form to indicate HST/WFC3 filters: WFC3/F098M $\rightarrow z_{\mathrm{098}}$; WFC3/F105W $\rightarrow Y_{\mathrm{105}}$; WFC3/F125W $\rightarrow J_{\mathrm{125}}$; WFC3/F140W $\rightarrow JH_{\mathrm{140}}$ and WFC3/F160W $\rightarrow H_{\mathrm{160}}$. We adopt a cosmology with $H_{0}=70$ km/s/Mpc, $\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 dataset 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 accompaining paper (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 (DR3111https://www.eso.org/sci/observing/phase3/data_releases/uvista_dr3.pdf) of the UltraVISTA program (McCracken et al. 2012). This data release provides mosaics in the $Y,J,H$ and $K_{\mathrm{S}}$ broad bands together with a narrow band centered at $1.18\mu$m (NB118). The mosaics in the broad band filters are characterized by four ultradeep stripes reaching $Y$-to-$K_{\mathrm{S}}\sim 24.8-25.3$ mag ($5\sigma,1\farcs 2$ aperture diameter corrected to total), alternating with four deep stripes ($Y$-to-$K_{\mathrm{S}}\sim 23.7-24.4$ mag, $5\sigma,1\farcs 2$ aperture to total), for a total area of $\sim 1.6$ degree2. 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 the Canada-France-Hawaii Legacy Survey (CFHTLS) and Subaru/Suprime-Cam (hereafter 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\mu$m and $4.5\mu$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 Hyper-Suprime-Cam (SPLASH, PI: Capak), and the complete set of observations of the Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS, PI: Caputi - Caputi et al. 2017; Ashby et al. 2017, in preparation).

Table 1 lists the $5\sigma$ depths of the adopted ground-based and *Spitzer/*IRAC mosaics. They were measured as the standard deviation of fluxes in $\sim 4000$ empty apertures of $1\farcs 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\mu$m and $4.5\mu$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 $\approx 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 $\sim 1.5$ lower than the standard deviation, thus corresponding to $\sim 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\sigma$ depths and recovered values within $0.1$ mag from those presented by Bowler et al. (2014)222This test was performed on $1\farcs 8$ empty apertures for full consistency with Bowler et al. (2014)..

Our search was carried out on the whole 1.6 degree2 of the UltraVISTA field. The mosaics of the optical and ground-based NIR bands were 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 SExtractor (Bertin & Arnouts 1996) in dual mode. Source detection was performed on the square root of the $\chi^{2}$ image (Szalay et al. 1999) built using the UltraVISTA $J,H$ and $K_{\mathrm{S}}$ band science and rms-map mosaics. Total fluxes were computed from $1\farcs 2$-diameter apertures and applying a correction based on the point-spread function (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 neighbouring objects within a radius of $15\farcs 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 neighbour-cleaned source. For this work, we adopted an aperture of $1\farcs 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 neighbouring sources.

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

[TABLE]

for the Lyman break, and

[TABLE]

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 $\geq 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 Eq. (1) includes two different Lyman break criteria: the first one selects galaxies whose Lyman break enters the Y band, i.e., whose redshift is $\gtrsim 7.5$ and the second one selects galaxies whose Lyman break enters the J band, i.e, when the redshift is $\gtrsim 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 vs 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\gtrsim 8$ galaxies brighter than $H=25.8$ mag.

Out of the 16 candidate galaxies at $z\gtrsim 8$, we selected five (labelled UVISTA-Y-1, UVISTA-Y-5, UVISTA-Y-6, UVISTA-J-1 and UVISTA-J-2) with plausible $z_{\mathrm{phot}}\gtrsim 8.5$ solutions, that stood out by their unprecedented brightness ($24.4\leq H\leq 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 March 27th, 2017 to March 29th, 2017. 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_{\mathrm{098}}$ hereafter), $\sim 460$ s (0.18 orbits) in the F125W band ($J_{\mathrm{125}}$ hereafter), and $\sim 410$ s (0.17 orbits) in the F160W band ($H_{\mathrm{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 2.

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 1 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\farcs 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 $\lesssim 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\farcs 2$ (diameter) and corrected to total using the light profile of each source. Remarkably, the optical data now include the mosaics from the HyperSuprimeCam Survey (HSC - Aihara et al. 2017a, b), not available at the time of the original selection of 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, which could introduce systematics in the estimates of the total fluxes. As a sanity check, we repeated the photometry with an aperture of $1\farcs 8$ diameter. The recovered total flux densities are on average within a few percent of the measurements based on the $1\farcs 2$ apertures, and within $\sim 15\%$ ($\sim 1\sigma$) in the worst cases. We therefore considered the photometry obtained adopting the $1\farcs 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 dataset.

For *HST/*WFC3 bands, we estimated the noise associated to the background from the dispersion of values in 200 $0\farcs 6$ apertures randomly placed across the image, free from sources according to the segmentation map, and repeating this process 20 times to increase the statistical significance. The final value was obtained 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 subtracting all the detected objects from the science cutout, was computed for apertures of $1\farcs 2$. As an additional step, the rms value was taken to be the maximum between the rms initially estimated by mophongo and the rms obtained from the empty apertures (see Sect. 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 uncertainties across different sources, in a given band, less immediate. Nonetheless, such method provides a robust and more comprehensive estimate of flux uncertainties.

One example for the above behaviour 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 $\sim 0.76$, which would correspond to an increase of the rms background for UVISTA-Y-6 by a factor of $\sim 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 ($\sim 2.5$ versus $\sim 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 neighbour subtraction. This, together with the additional information provided by the HSC survey (especially for UVISTA-Y-1 which lacks coverage from the CFHTLS survey), 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\sigma$ 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_{\mathrm{V}}=2.0,5.0,8.0$ mag curves.

The input catalog consisted in the flux measurements listed in Table 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_{\mathrm{phot}}=8.38^{+0.35}_{-0.43},8.74^{+0.45}_{-0.47}$ and $8.53^{+0.53}_{-0.80}$, 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_{\mathrm{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$ (grey 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\sim 8$ best-fit template. Specifically, the brown dwarf template is inconsistent with observations in the IRAC $3.6\mu$m and $4.5\mu$m bands, and, for UVISTA-Y-1 also with the $JH_{\mathrm{140}}$ measurement. This is reflected by the poorer best-fit $\chi^{2}$’s, with $\chi^{2}=54.8,58.5$ and $41.2$, respectively, for the 3 sources. Forcing the solution to be at $z<6$ results in $z_{\mathrm{phot}}\sim 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 solution. However, it remains in tension with the data, resulting in $\chi^{2}$’s 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 template for UVISTA-J-1 and UVISTA-J-2, displayed in Figure 4, are inconsistent with our measured fluxes in the IRAC $3.6\mu$m and $4.5\mu$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 to 1.5 randomly extracted from a uniform distribution. The new catalog was analysed 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 Sect. 2 (typical differences were factors of 1.5). We found photometric redshifts $z_{\mathrm{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$ 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 mis-interpreted 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\mu$m, $4.5\mu$m and $5.8\mu$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\mu$m data as it is likely not contaminated by strong nebular emission yet it provides constraints for the SED modelling. Indeed, the $K_{\mathrm{S}}-[4.5]$ color could be interpreted as the Balmer break, guiding the fit towards higher redshift solutions. These new measurements resulted in lower photometric redshifts: $z_{\mathrm{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_{\mathrm{phot}}=10.1^{+1.4}_{-0.8}$, consistent with the initial selection. The new $H_{160}$ observations point to a much redder overall near-IR 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_{\mathrm{phot}}=2.2^{+0.6}_{-0.4}$). After further inspection, we conclude that the likely reason for this is a higher flux measurement in the $3.6\mu$m band we obtained using the new HST dataset as morphological and positional prior which allowed for a more accurate subtraction of the neighbours, 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-Rays catalogs from XMM and Chandra (Cappelluti et al. 2009 and Civano et al. 2016; Marchesi et al. 2016, respectively). Finally, visual inspection of the MIPS $\mu$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 $\sim 0\farcs 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\sigma$) in the HST/WFC3 $z_{\mathrm{098}}$ band, strengthening the evidence that this is a $z>7.5$ LBG. The WFC3 photometry in the $J_{\mathrm{125}}$ and $H_{\mathrm{160}}$ is consistent at $1\sigma$ or better with that in the UltraVISTA $J$ and $H$ bands, respectively. This source also benefits from additional WFC3 coverage in the $Y_{\mathrm{105}}$ and $JH_{\mathrm{140}}$ bands from the SUSHI program (PI: N. Suzuki). The measurement in the $JH_{\mathrm{140}}$ band is consistent with the best-fit SED. The $p(z)$ is characterized by a solution at $8.02^{+0.41}_{-0.49}$, with a marginal secondary peak at $z\sim 1.8$ ($p(z<6)=0.12$). The $Y_{\mathrm{105}}$ band shows a $2.2\sigma$ detection, as expected if the source is at $z\sim 8$.

UVISTA-Y-5: The source is undetected ($<2\sigma$) in the HST/WFC3 $z_{\mathrm{098}}$ band, strongly favouring a $z>7.5$ solution for this source. The WFC3 photometry in the $J_{\mathrm{125}}$ is consistent at $1\sigma$ with that in the UltraVISTA $J$ band. However, the flux measurement in the $H_{\mathrm{160}}$ results in a $1.42\sigma$ detection only, and is consistent with the UltraVISTA H-band photometry at $\sim 3\sigma$ level. In Appendix A we analyse 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 at $z\sim 1.8$ ($p(z<6)=0.009$).

UVISTA-Y-6: This source is undetected ($<2\sigma$) in the HST/WFC3 $z_{\mathrm{098}}$ band, again favoring a $z>7.5$ solution. The WFC3 photometry in the $J_{\mathrm{125}}$ and $H_{\mathrm{160}}$ is consistent at $1\sigma$ or better with that in the UltraVISTA $J$ and $H$ bands, respectively. The $p(z)$ shows a distribution with best fit solution $z_{\mathrm{phot}}=8.35^{+0.65}_{-0.81}$ 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_{\mathrm{160}}$ band ($<2\sigma$), making it consistent with the UltraVISTA $H$ band only at $\sim 3.5\sigma$ 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_{\mathrm{125}}$ band, instead, is consistent with that in the UltraVISTA $J$ band at $1\sigma$ level. The lower-z solution is enforced by the fact that the source is detected in the HSC $z$ band at $1.7\sigma$ level. The fiducial photometric redshift is $2.05^{+0.49}_{-0.46}$. 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_{\mathrm{125}}$ and $H_{\mathrm{160}}$ is consistent with that in the corresponding UltraVISTA bands at $\sim 1.8-3\sigma$, respectively, with nominal redshift of $1.96^{+7.04}_{-0.33}$. Similarly to what done for UVISTA-Y-5 and UVISTA-J-1, in Appendix A we analyse 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_{\mathrm{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 the present. 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 $\sim 1.8$ mag more luminous (a factor $\sim 5.3\times$) than the estimated characteristic magnitude at $z\sim 9$.

5.2 $\beta-M_{\mathrm{UV}}$ relation

We measured the rest-frame UV slope ($\beta$) fitting a power-law of the form $f_{\lambda}\propto\lambda^{\beta}$ to the fluxes in the $H_{\mathrm{160}}$ and in the UltraVISTA $H$ and $K_{\mathrm{S}}$ bands. The results are presented in Figure 6 and listed in Table 4. These slopes have an average value of $\beta=-1.91\pm 0.26$ and are consistent with the recent determination of the UV slopes of Bowler et al. (2017) for LBGs with similar luminosity ($M_{\mathrm{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_{\mathrm{UV}}\sim-21$) $z\sim 10$ LBGs by Wilkins et al. (2016). For comparison, in the plot we also show the bi-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; Finkelstein et al. 2012; Castellano et al. 2012; Dunlop et al. 2013; Jiang et al. 2013; Bouwens et al. 2014; Duncan et al. 2014; Rogers et al. 2014; Duncan & Conselice 2015; see also Oesch et al. 2013 and Stefanon et al. 2016 for similar relations of $\beta$ and rest-frame optical luminosities). This behaviour 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 $\beta-M_{\mathrm{UV}}$ 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 $\beta$ measurements make them consistent at $<1\sigma$ with those relations, thus preventing from further inspecting any differential evolutionary path of $\beta$ 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 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 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 signal-to-noise in most of the WFC3 images, for our analysis we only considered the $JH_{\mathrm{140}}$ band of UVISTA-Y-1 ($\sim$ rest-frame 1600 Å), i.e., the highest signal-to-noise 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 signal to noise 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 and Bowler et al. 2017). We then verified that $R_{e}$ does not systematically change ($\leq 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 neighbours within $5\farcs 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,\mathrm{circ}}=R_{e}\sqrt{b/a}$, where $b/a$ is the minor-to-major axis ratio. As a consistency check, we also derive $R_{e}$ using SExtractor. In this case, the final value for $R_{e,\mathrm{SE}}=\sqrt{R_{e,\mathrm{obs,SE}}^{2}-r_{\mathrm{PSF}}^{2}}$, with $R_{e,\mathrm{obs,SE}}$ the effective radius measured by SExtractor and $R_{\mathrm{PSF}}$ that of the $JH_{140}$ PSF, with $R_{\mathrm{PSF}}=0\farcs 12$.

We find $R_{e,\mathrm{circ}}=0.9\pm 0.2$ kpc from galfit, consistent with $R_{e,\mathrm{SE}}=0.7$ kpc estimated with SExtractor. In Table 5 we list the main morphological parameter we obtain from the two methods. Our values are consistent at $1\sigma$ with estimates of $R_{e}$ for $M_{\mathrm{UV}}\sim-22$ LBGs at $z\sim 7$ ($R_{e,\mathrm{circ}}=0.6-0.9\pm 0.2$ from a stacking analysis - Bowler et al. 2017) and $z\sim 11$ ($R_{e,\mathrm{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_{\mathrm{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 star-formation rate (SFR) surface density ($\Sigma_{\mathrm{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,\mathrm{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_{\mathrm{SFR}}\sim 11M_{\odot}$ yr*-1* 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_{\mathrm{SFR}}$ for $L\lesssim L^{*}$ galaxies a factor $\sim 3\times$ lower).

Interestingly, but unsurprisingly, $\Sigma_{\mathrm{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_{\mathrm{SFR}}-z$ relation (e.g., Oesch et al. 2010; Ono et al. 2013). This is qualitatively consistent with the increase of the star-formation rate density with cosmic time combined with the increase in size. Our (mildly) evolving $\Sigma_{\mathrm{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 this could be the case at even higher redshift, data with better S/N is 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 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 include the full suite of HST, ground-based, and Spitzer/IRAC images. For the ground-based and Spitzer/IRAC data the mock sources were convolved with appropriate kernels to match the lower resolution PSF. To simulate IRAC colors we assume a continuum flat in $f_{\nu}$ and emission lines with fixed rest-frame EW($H\alpha$+[NII]+[SII]) = 300Å and rest-frame EW([OIII]+$H\beta$) = 500Å, consistent with the results of Labbé et al. (2013); Stark et al. (2013); Smit et al. (2014, 2015); Rasappu et al. (2016). The same detection and selection criteria as described in Sect. 2 were then applied to the simulated images to recover the completeness as a function of magnitude and the selection 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 ultra-deep and deep stripes, used to estimate the co-moving 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 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_{\mathrm{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_{\mathrm{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_{\mathrm{max}}$ method directly provides the normalisation 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 obtain a volume density of $\Phi=8.49^{+8.23}_{-4.60}\times 10^{-7}$ Mpc*-3* mag*-1* at $M_{\mathrm{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, excluding UV LFs from samples based on galaxy cluster (Bradley et al. 2012; McLure et al. 2013; Schenker et al. 2013; Schmidt et al. 2014; Bouwens et al. 2015; Finkelstein et al. 2015). Our measurement constitutes the first volume density estimate for $M_{\mathrm{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 Labbé et al. (2017, in preparation) includes three more potential galaxies at $7.5\lesssim z_{\mathrm{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_{\mathrm{UV}}=-22.21\pm 0.25$ the volume density is $\Phi=17.3^{+10.3}_{-6.9}\times 10^{-7}$ Mpc*-3* 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 parameterised 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 Sect. 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, though, as the corresponding absolute magnitude lies at the intersection of the Schechter and the DPL forms.

Because the $z_{\mathrm{phot}}$ solutions for half of our sample of candidate $z\gtrsim 8$ LBGs may have values close to $z_{\mathrm{phot}}\sim 9$ when including the IRAC bands, we also computed the volume density associated to the three sources (UVISTA-Y-5, UVISTA-Y-6 and UVISTA-Y-2) with $z_{\mathrm{phot}}>8.5$. We obtain $\Phi=8.5^{+8.2}_{-4.6}\times 10^{-7}$ Mpc*-3* 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. (2016b); Calvi et al. (2016); McLeod et al. (2016); 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. (2016b). In the same panel we also plot the bright-end of the dual power-law 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 the 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 the volume density estimate 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 hundreds 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 SFR.

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 Labbé et al. (2017, in 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 (labelled 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_{\mathrm{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 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 Labbé et al. (2017, in preparation).

The HST/WFC3 observations were performed in the $z_{\mathrm{098}}$, $J_{\mathrm{125}}$ and $H_{\mathrm{160}}$ bands (Figure 1) in March 2017, 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_{\mathrm{105}}$ and $JH_{\mathrm{140}}$ bands to $25.6$ and $25.5$ mag, respectively ($5\sigma$, $1\farcs 2$ aperture diameter - Figure 2). 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, b).

Our analysis confirms the photometric redshift of three sources (UVISTA-Y-1, UVISTA-Y-5 and UVISTA-Y-6) to be $8.0<z_{\mathrm{photo}}<8.7$ (Figure 3). Their measured luminosity $M_{\mathrm{UV}}\sim-22.3$ makes them perhaps the brightest, most reliable galaxies at $z\sim 8-9$ identified to date. The uniquely deep optical, near-IR, 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_{\mathrm{phot}}\sim 2$ (Figure 4).

The three $z\sim 8$ candidate LBGs are characterized by average UV continuum slopes $\beta=-1.91\pm 0.26$, 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 from deriving any further conclusion on differential evolution.

For our bright source UVISTA-Y-1 we measure a size of $R_{e,\mathrm{circ}}=0.9\pm 0.2$ 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_{\mathrm{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^{+8.23}_{-4.60}\times 10^{-7}$ Mpc*-3* mag*-1* at $M_{\mathrm{UV}}=-22.21\pm 0.25$. This constitutes the first measurement of the number density of $M_{\mathrm{UV}}<-22$ mag LBGs at $z\sim 8$ based on actual detection of sources. We also estimate the volume density of star-forming galaxies at $z\sim 8$, including the full sample of six $M_{\mathrm{UV}}\lesssim-22$ mag galaxies identified in 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_{\mathrm{UV}}\lesssim-22$ mag. The measured volume density associated to this complete sample is $\Phi=17.3^{+10.3}_{-6.9}\times 10^{-7}$ Mpc*-3* mag*-1* at $M_{\mathrm{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.

This work further stresses the importance of the high-resolution imaging provided by HST in the study of the galaxy populations at early cosmic epochs, enabling to refine the photometric redshifts and indentify interlopers, resulting in cleaner samples. Our results, however, are based on photometric redshifts from broad-band 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 analysed in this work offer two possibilities: 1) they constitute prime candidates for future spectroscopic follow-up with JWST; 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 operations of JWST.

Acknowledgements: MS would like to thank Pieter van Dokkum, Jorryt Matthee and Corentin Schreiber for useful discussion. KIC acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. 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 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, Queen s 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 Science de l’Univers of the Centre National 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 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.

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 favouring 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\farcs 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 zeropoint offset is already included in the photometry we provide in Table 3.

In principle, the different releases of the UltraVISTA dataset 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 CIII]$\lambda\lambda 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 CIII] at $z\sim 7$). Even assuming that CIII] is dominating the flux in the UltraVISTA $H$ band, similar EW for CIV would be required to explain the flux in the $K_{\mathrm{S}}$ band, together with a very red SED red-ward of $K_{\mathrm{S}}$ to match the fluxes in the $3.6\mu$m and $4.5\mu$m bands. One further possibility could be that the observed discrepancy was originated by [OIII]+H$\beta$ nebular lines if the galaxy had a redshift $z\sim 2.4-2.6$. In this case, EW0([OIII]+H$\beta$)$\gtrsim 1300$ Å. This value is larger than inferred from conversions of observed H$\alpha$ 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 ) neighbouring star located $\sim 40$ arcsec north-east. Even though the procedure we adopted to extract the flux measurements takes care of estimating the background, we can not 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\geq 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. 2016a).

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 Labbé et al. (2017, in preparation), that we included in our LF estimates (Sect. 5.4). Their selection followed the same methods described in Sect. 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., neighbouring, 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 of detecting the source of interest.

In Figure 10 we present image cutouts in the optical, NIR and IRAC $3.6\mu$m and $4.5\mu$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 $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_{\mathrm{UV}}\lesssim-22$ mag.