The Brightest $z\gtrsim8$ Galaxies over the COSMOS UltraVISTA Field

TL;DR

This paper reports the discovery of 16 ultrabright galaxy candidates at redshift around 8 in the COSMOS/UltraVISTA field, using deep multi-wavelength data and strict selection criteria, and analyzes their properties and implications for the UV luminosity function.

Contribution

First identification of 16 ultrabright z~8 galaxy candidates using ground-based and Spitzer data, with detailed analysis of their properties and the bright end of the UV luminosity function.

Findings

Sample includes 18 galaxies, 10 with high probability of being at z>7.

Bright z~8 galaxies have bluer UV slopes than z~6 and z~7 counterparts.

The UV luminosity function at z~8 can be modeled by Schechter or double power-law forms.

Abstract

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

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 three 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.2$	$26.7$
SSC $B$	$1.7$	$27.4$
HSC $g$ccThe HyperSuprimeCam 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.	$2.1$	$26.7$
CFHTLS $g$	$2.1$	$26.8$
SSC $V$	$2.1$	$26.4$
HSC $r$ccThe HyperSuprimeCam 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.	$1.7$	$26.8$
CFHTLS $r$	$2.0$	$26.4$
SSC $r^{+}$	$2.0$	$26.6$
SSC $i^{+}$	$1.9$	$26.2$
CFHTLS $y$	$1.9$	$26.1$
CFHTLS $i$	$1.9$	$26.0$
HSC $i$ccThe HyperSuprimeCam 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.	$1.8$	$26.3$
CFHTLS $z$	$2.0$	$25.2$
HSC $z$ccThe HyperSuprimeCam 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.	$1.7$	$25.9$
SSC $z^{+}$	$2.2$	$25.0$
HSC $y$ccThe HyperSuprimeCam 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.	$2.1$	$24.9$
UVISTA $Y$	$2.5$	$25.4/24.5$
UVISTA $J$	$2.3$	$25.4/24.4$
UVISTA $H$	$2.2$	$25.1/24.1$
UVISTA $K_S$	$2.1$	$24.8/23.7$
IRAC $3.6\mu$m	$2.7\text{<span><sup>d</sup><span><span><sup>d</sup>Aperture corrections for IRAC bands refer to the <xmath><xmapp><xmtok>⁢</xmtok><xmtok>1</xmtok><xmdual><xmtok><xmapp><xmtok><xmtok>.</xmtok><xmwrap><xmtok>′</xmtok><xmtok>′</xmtok></xmwrap></xmtok></xmapp></xmtok></xmdual><xmtok>8</xmtok></xmapp></xmath> diameter.</span></span></span>}$	$25.4/24.9/24.5$
IRAC $4.5\mu$m	$2.7\text{<span><sup>d</sup><span><span><sup>d</sup>Aperture corrections for IRAC bands refer to the <xmath><xmapp><xmtok>⁢</xmtok><xmtok>1</xmtok><xmdual><xmtok><xmapp><xmtok><xmtok>.</xmtok><xmwrap><xmtok>′</xmtok><xmtok>′</xmtok></xmwrap></xmtok></xmapp></xmtok></xmdual><xmtok>8</xmtok></xmapp></xmath> diameter.</span></span></span>}$	$25.3/24.7/24.3$
IRAC $5.8\mu$m	$3.4\text{<span><sup>d</sup><span><span><sup>d</sup>Aperture corrections for IRAC bands refer to the <xmath><xmapp><xmtok>⁢</xmtok><xmtok>1</xmtok><xmdual><xmtok><xmapp><xmtok><xmtok>.</xmtok><xmwrap><xmtok>′</xmtok><xmtok>′</xmtok></xmwrap></xmtok></xmapp></xmtok></xmdual><xmtok>8</xmtok></xmapp></xmath> diameter.</span></span></span>}$	$20.8$
IRAC $8.0\mu$m	$4.1\text{<span><sup>d</sup><span><span><sup>d</sup>Aperture corrections for IRAC bands refer to the <xmath><xmapp><xmtok>⁢</xmtok><xmtok>1</xmtok><xmdual><xmtok><xmapp><xmtok><xmtok>.</xmtok><xmwrap><xmtok>′</xmtok><xmtok>′</xmtok></xmwrap></xmtok></xmapp></xmtok></xmdual><xmtok>8</xmtok></xmapp></xmath> diameter.</span></span></span>}$	$20.6$

Table 2. Table 2: Candidate z ∼ 8 similar-to 𝑧 8 z\sim 8 LBGs with HST /WFC3 F160W coverage

ID	PID	PI	Depth
			[mag]
UVISTA-Y1	14895	R. Bouwens	$24.7$
UVISTA-Y2	14114	P. van Dokkum	$24.9$
UVISTA-Y3aaHST/WFC3 imaging suggests this source is potentially multiple. See Sect. 5.1 for details.	13868	D. Kocevski	$26.5$
UVISTA-Y4	14114	P. van Dokkum	$24.9$
UVISTA-Y5	14895	R. Bouwens	$24.9$
UVISTA-Y6	14895	R. Bouwens	$25.0$
UVISTA-Y7	14114	P. van Dokkum	$24.9$
UVISTA-Y8	13641	P. Capak	$25.7$
UVISTA-Y9	14114	P. van Dokkum	$24.8$
UVISTA-Y10	14114	P. van Dokkum	$24.7$
UVISTA-Y11	12440	S. Faber	$26.6$
UVISTA-Y13	14114	P. van Dokkum	$24.9$
UVISTA-Y14	14114	P. van Dokkum	$24.8$
UVISTA-Y16	14114	P. van Dokkum	$24.7$

Table 3. Table 3: Sample of candidate z ∼ 8 similar-to 𝑧 8 z\sim 8 LBGs

ID	R.A.	Dec.	$m_H$aa$H$-band magnitude and associated $1\sigma$ uncertainty estimated from the UltraVISTA DR3 mosaic.	$Y-J$bbUpper/lower limits to be intended as $1\sigma$. In computing these colors, we replaced negative fluxes with their corresponding $1\sigma$ uncertainty. See Tables 7, 8 and 9 for a complete listing of flux densities in all bands.	$\left[3.6\right]-\left[4.5\right]$bbUpper/lower limits to be intended as $1\sigma$. In computing these colors, we replaced negative fluxes with their corresponding $1\sigma$ uncertainty. See Tables 7, 8 and 9 for a complete listing of flux densities in all bands.	$z_{\mathrm{phot}}$ccPhotometric redshift and 68% confidence interval of the best-fitting template from EAzY.
	[J2000]	[J2000]	[mag]	[mag]	[mag]
UVISTA-Y1d,*d,*footnotemark:	$09:57:47.900$	$+02:20:43.66$	$24.8\pm 0.1$	$>2.1$	$0.4\pm 0.2$	${8.53}_{-0.62}^{+0.53}$
UVISTA-Y2**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:02:12.558$	$+02:30:45.71$	$24.8\pm 0.2$	$>2.2$	$0.5\pm 0.1$	${8.21}_{-0.49}^{+0.50}$
UVISTA-Y3aeeThese candidate LBGs were initially identified as a single source on the UltraVISTA NIR bands. Successive analysis including COSMOS/DASH suggests these are three distinct objects. The corresponding observables when a single object is assumed are: R.A.= 10:00:32.322; Dec=1:44:31.26, $m_H=25.0\pm 0.1$ mag; $Y-J=1.1\pm 0.4$ mag; $\left[3.6\right]-\left[4.5\right]=0.3\pm 0.1$ mag and $z_{\mathrm{phot}}={7.62}_{-0.28}^{+0.14}$	$10:00:32.324$	$+01:44:30.86$	$25.5\pm 0.3$	$>0.9$	$0.6\pm 0.5$	${8.68}_{-1.21}^{+0.93}$
UVISTA-Y3beeThese candidate LBGs were initially identified as a single source on the UltraVISTA NIR bands. Successive analysis including COSMOS/DASH suggests these are three distinct objects. The corresponding observables when a single object is assumed are: R.A.= 10:00:32.322; Dec=1:44:31.26, $m_H=25.0\pm 0.1$ mag; $Y-J=1.1\pm 0.4$ mag; $\left[3.6\right]-\left[4.5\right]=0.3\pm 0.1$ mag and $z_{\mathrm{phot}}={7.62}_{-0.28}^{+0.14}$	$10:00:32.317$	$+01:44:31.48$	$26.1\pm 0.5$	$>0.9$	ffThis IRAC color is based on flux density estimate in both bands.	${8.90}_{-1.18}^{+1.24}$
UVISTA-Y3ceeThese candidate LBGs were initially identified as a single source on the UltraVISTA NIR bands. Successive analysis including COSMOS/DASH suggests these are three distinct objects. The corresponding observables when a single object is assumed are: R.A.= 10:00:32.322; Dec=1:44:31.26, $m_H=25.0\pm 0.1$ mag; $Y-J=1.1\pm 0.4$ mag; $\left[3.6\right]-\left[4.5\right]=0.3\pm 0.1$ mag and $z_{\mathrm{phot}}={7.62}_{-0.28}^{+0.14}$	$10:00:32.350$	$+01:44:31.73$	$26.0\pm 0.5$	$>-0.5$	$0.7\pm 0.5$	${9.29}_{-2.10}^{+1.58}$
UVISTA-Y4**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:00:58.485$	$+01:49:55.96$	$24.9\pm 0.2$	$1.0\pm 0.4$	$0.1\pm 0.2$	${7.42}_{-0.20}^{+0.19}$
UVISTA-Y5d,*d,*footnotemark:	$10:00:31.886$	$+01:57:50.23$	$24.9\pm 0.2$	$>1.3$	$0.8\pm 0.3$	${8.60}_{-0.65}^{+0.58}$
UVISTA-Y6ddThese sources were already presented in Stefanon et al. (2017b). We propose them here again for completeness, noting that their associated parameters in the present work were computed excluding the information from the HST bands. We refer the reader to Stefanon et al. (2017b) for a more complete analysis.	$10:00:12.506$	$+02:03:00.50$	$25.3\pm 0.3$	$>1.5$	$0.3\pm 0.4$	${8.32}_{-0.92}^{+0.66}$
UVISTA-Y7**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$09:59:02.566$	$+02:38:06.05$	$25.5\pm 0.4$	$>1.3$	$\cdots$$\dagger$$\dagger$footnotemark:	${8.47}_{-0.73}^{+0.72}$
UVISTA-Y8**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:00:47.544$	$+02:34:04.84$	$25.4\pm 0.3$	$>1.4$	$1.0\pm 0.8$	${8.34}_{-0.58}^{+0.60}$
UVISTA-Y9	$09:59:09.621$	$+02:45:09.68$	$25.4\pm 0.3$	$0.8\pm 0.7$	$\cdots$$\dagger$$\dagger$footnotemark:	${7.69}_{-0.71}^{+0.99}$
UVISTA-Y10**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:01:47.495$	$+02:10:15.37$	$25.3\pm 0.3$	$>1.6$	$0.9\pm 0.7$	${8.25}_{-0.60}^{+0.61}$
UVISTA-Y11**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:00:19.607$	$+02:14:13.15$	$25.2\pm 0.3$	$>1.4$	$0.8\pm 0.4$	${8.64}_{-0.72}^{+0.66}$
UVISTA-Y12**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:00:15.975$	$+02:43:32.96$	$25.6\pm 0.4$	$>1.2$	$0.2\pm 0.8$ffThis IRAC color is based on flux density estimate in both bands.	${8.70}_{-0.74}^{+0.61}$
UVISTA-Y13	$09:58:45.561$	$+01:53:41.79$	$25.8\pm 0.4$	$>1.1$	$0.8\pm 0.7$	${8.54}_{-1.18}^{+0.79}$
UVISTA-Y14	$10:00:12.568$	$+01:54:28.50$	$25.6\pm 0.4$	$>1.1$	$0.1\pm 0.6$	${7.55}_{-2.68}^{+1.71}$
UVISTA-Y15	$09:57:35.795$	$+02:11:57.81$	$25.6\pm 0.4$	$1.1\pm 0.9$	f,gf,gfootnotemark:	${7.64}_{-1.13}^{+1.13}$
UVISTA-Y16**These sources have a probability $p\left(z>7\right)\ge 0.97$, suggesting these may be fairly robust candidates of bright LBGs.	$10:01:56.333$	$+02:34:16.25$	$25.3\pm 0.3$	$1.2\pm 0.7$	$0.6\pm 0.4$	${7.90}_{-0.57}^{+0.74}$

Table 4. Table 4: Main physical parameters for the sample of candidate z ∼ 8 similar-to 𝑧 8 z\sim 8 LBGs

ID	$M_{\mathrm{UV}}$	UV slope $\beta$	$u-g$	$\log \left(M_{\star }\right)$	$\log \left(\mathrm{SFR}\right)$	$\log \left(\mathrm{sSFR}\right)$	$\log \left(\mathrm{age}\right)$	$A_V$
	[mag]		[mag]	[$M_{\odot }$]	[$M_{\odot }{\mathrm{yr}}^{-1}$]	[yr-1]	[yr]	[mag]
UVISTA-Y1	$-22.48\pm 0.15$	$-{1.5}_{-0.7}^{+0.4}$	$0.45\pm 0.18$	${10.0}_{-0.4}^{+0.9}$	${1.59}_{-9.55}^{+1.02}$	$-{8.4}_{-9.8}^{+1.8}$	${7.30}_{-0.61}^{+1.42}$	${0.9}_{-0.9}^{+0.0}$
UVISTA-Y2	$-22.37\pm 0.20$	$-{2.6}_{-0.5}^{+0.5}$	$0.66\pm 0.21$	${9.0}_{-1.2}^{+0.3}$	${1.98}_{-7.57}^{+0.65}$	$-{7.0}_{-8.6}^{+0.8}$	${7.00}_{-0.50}^{+1.80}$	${0.4}_{-0.4}^{+0.3}$
UVISTA-Y3aaaThese three candidate LBGs were originally identified as a single source, successively de-blended using data from the COSMOS/DASH program (see Sect. 5.1 and Figure 10). When we do not deblend the source, we obtain $M_{\mathrm{UV}}=-22.00\pm 0.16$ mag, $\beta =-1.8\pm 0.7$, $u-g=0.58\pm 0.16$ mag, $\log \left(M_{\star }/M_{\odot }\right)={9.9}_{-0.3}^{+0.6}$, $\log \left(\mathrm{SFR}/M_{\odot }/{\mathrm{yr}}^{-1}\right)={1.63}_{-3.77}^{+0.38}$, $\log \left(\mathrm{sSFR}/{\mathrm{yr}}^{-1}\right)=-{8.2}_{-3.8}^{+0.9}$, $\log \left(\mathrm{age}/\mathrm{yr}\right)={8.20}_{-1.16}^{+0.60}$ and $A_V={0.5}_{-0.5}^{+0.5}$ mag.	$-21.77\pm 0.32$	$-{1.5}_{-0.9}^{+0.7}$	$0.57\pm 0.39$	${9.8}_{-0.3}^{+1.3}$	$-{1.34}_{-7.01}^{+4.06}$	$-{11.1}_{-7.0}^{+4.9}$	${8.00}_{-1.50}^{+0.80}$	${0.0}_{-0.0}^{+1.1}$
UVISTA-Y3baaThese three candidate LBGs were originally identified as a single source, successively de-blended using data from the COSMOS/DASH program (see Sect. 5.1 and Figure 10). When we do not deblend the source, we obtain $M_{\mathrm{UV}}=-22.00\pm 0.16$ mag, $\beta =-1.8\pm 0.7$, $u-g=0.58\pm 0.16$ mag, $\log \left(M_{\star }/M_{\odot }\right)={9.9}_{-0.3}^{+0.6}$, $\log \left(\mathrm{SFR}/M_{\odot }/{\mathrm{yr}}^{-1}\right)={1.63}_{-3.77}^{+0.38}$, $\log \left(\mathrm{sSFR}/{\mathrm{yr}}^{-1}\right)=-{8.2}_{-3.8}^{+0.9}$, $\log \left(\mathrm{age}/\mathrm{yr}\right)={8.20}_{-1.16}^{+0.60}$ and $A_V={0.5}_{-0.5}^{+0.5}$ mag.	$-21.23\pm 0.54$	$-{3.2}_{-0.0}^{+2.2}$	$-0.22\pm 1.82$	${8.7}_{-0.0}^{+0.1}$	$-{0.28}_{-0.07}^{+0.00}$	$-{9.0}_{-0.0}^{+0.0}$	${7.40}_{-0.02}^{+0.00}$	${0.0}_{-0.0}^{+0.0}$
UVISTA-Y3caaThese three candidate LBGs were originally identified as a single source, successively de-blended using data from the COSMOS/DASH program (see Sect. 5.1 and Figure 10). When we do not deblend the source, we obtain $M_{\mathrm{UV}}=-22.00\pm 0.16$ mag, $\beta =-1.8\pm 0.7$, $u-g=0.58\pm 0.16$ mag, $\log \left(M_{\star }/M_{\odot }\right)={9.9}_{-0.3}^{+0.6}$, $\log \left(\mathrm{SFR}/M_{\odot }/{\mathrm{yr}}^{-1}\right)={1.63}_{-3.77}^{+0.38}$, $\log \left(\mathrm{sSFR}/{\mathrm{yr}}^{-1}\right)=-{8.2}_{-3.8}^{+0.9}$, $\log \left(\mathrm{age}/\mathrm{yr}\right)={8.20}_{-1.16}^{+0.60}$ and $A_V={0.5}_{-0.5}^{+0.5}$ mag.	$-21.37\pm 0.53$	$-{2.0}_{-0.0}^{+1.7}$	$0.81\pm 0.57$	${10.2}_{-0.5}^{+1.6}$	${1.69}_{-28.23}^{+1.86}$	$-{8.5}_{-28.1}^{+2.3}$	${8.70}_{-2.20}^{+0.10}$	${0.8}_{-0.8}^{+1.1}$
UVISTA-Y4	$-22.11\pm 0.24$	$-{2.7}_{-0.4}^{+0.7}$	$0.47\pm 0.26$	${9.9}_{-0.2}^{+0.5}$	${1.23}_{-6.97}^{+0.61}$	$-{8.6}_{-7.0}^{+0.9}$	${8.50}_{-1.31}^{+0.30}$	${0.0}_{-0.0}^{+0.7}$
UVISTA-Y5	$-22.34\pm 0.24$	$-{1.7}_{-0.9}^{+0.8}$	$0.63\pm 0.27$	${9.0}_{-1.1}^{+0.4}$	${1.99}_{-7.68}^{+0.55}$	$-{7.0}_{-8.6}^{+0.8}$	${7.30}_{-0.80}^{+1.40}$	${0.4}_{-0.4}^{+0.2}$
UVISTA-Y6	$-21.92\pm 0.26$	$-{1.7}_{-0.8}^{+0.7}$	$0.49\pm 0.32$	${9.7}_{-0.5}^{+1.1}$	${1.36}_{-12.70}^{+1.33}$	$-{8.4}_{-12.9}^{+2.2}$	${7.30}_{-0.80}^{+1.50}$	${0.9}_{-0.9}^{+0.3}$
UVISTA-Y7	$-21.74\pm 0.36$	$-{2.0}_{-0.5}^{+0.7}$	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:
UVISTA-Y8	$-21.76\pm 0.35$	$-{2.8}_{-0.4}^{+0.9}$	$0.91\pm 0.45$	${8.3}_{-1.4}^{+0.1}$	${1.90}_{-1.41}^{+0.35}$	$-{6.4}_{-2.6}^{+0.2}$	${6.50}_{-0.00}^{+2.29}$	${0.0}_{-0.0}^{+0.5}$
UVISTA-Y9	$-21.66\pm 0.34$	$-{2.6}_{-0.6}^{+0.9}$	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:	$\cdots$$\dagger$$\dagger$footnotemark:
UVISTA-Y10	$-21.89\pm 0.31$	$-{2.2}_{-0.7}^{+1.1}$	$0.56\pm 0.39$	${8.3}_{-1.4}^{+0.0}$	${1.80}_{-4.34}^{+0.47}$	$-{6.5}_{-5.6}^{+0.3}$	${6.70}_{-0.20}^{+2.10}$	${0.0}_{-0.0}^{+0.5}$
UVISTA-Y11	$-22.04\pm 0.26$	$-{1.8}_{-1.3}^{+0.5}$	$0.67\pm 0.30$	${8.7}_{-1.2}^{+0.4}$	${1.76}_{-7.62}^{+0.60}$	$-{7.0}_{-8.6}^{+0.8}$	${7.30}_{-0.80}^{+1.44}$	${0.3}_{-0.3}^{+0.2}$
UVISTA-Y12	$-21.66\pm 0.40$	$-{2.1}_{-0.8}^{+1.3}$	$0.22\pm 0.64$	${9.1}_{-0.4}^{+0.9}$	${0.17}_{-2.88}^{+2.22}$	$-{9.0}_{-3.1}^{+2.8}$	${7.40}_{-0.90}^{+1.30}$	${0.2}_{-0.2}^{+0.3}$
UVISTA-Y13	$-21.39\pm 0.42$	$-{1.0}_{-0.7}^{+0.7}$	$0.50\pm 0.51$	${9.8}_{-0.3}^{+1.3}$	${0.70}_{-9.08}^{+1.82}$	$-{9.1}_{-9.1}^{+2.8}$	${7.50}_{-0.96}^{+1.28}$	${0.8}_{-0.8}^{+0.3}$
UVISTA-Y14	$-21.44\pm 0.40$	$-{3.0}_{-0.1}^{+1.7}$	$0.63\pm 0.54$	${9.3}_{-0.4}^{+1.2}$	${0.52}_{-9.14}^{+2.03}$	$-{8.8}_{-9.3}^{+2.6}$	${8.20}_{-1.70}^{+0.63}$	${0.0}_{-0.0}^{+0.8}$
UVISTA-Y15	$-21.50\pm 0.37$	$-{2.7}_{-0.4}^{+0.9}$	$-0.04\pm 0.68$	${8.8}_{-0.0}^{+0.2}$	$-{0.16}_{-0.02}^{+0.41}$	$-{9.0}_{-0.0}^{+0.6}$	${7.40}_{-0.10}^{+0.01}$	${0.0}_{-0.0}^{+0.0}$
UVISTA-Y16	$-21.80\pm 0.33$	$-{2.2}_{-0.8}^{+0.8}$	$0.39\pm 0.40$	${8.6}_{-0.9}^{+0.3}$	${1.62}_{-0.98}^{+0.56}$	$-{7.0}_{-1.8}^{+0.8}$	${7.30}_{-0.80}^{+1.50}$	${0.1}_{-0.1}^{+0.4}$

Table 5. Table 5: Observed and rest-frame properties for candidate z ∼ 8 similar-to 𝑧 8 z\sim 8 galaxies identified in the UltraVISTA DR3 observations

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

Table 6. Table 6: V max determinations of the UV LF

$M_{\mathrm{UV}}$	$\varphi$
[mag]	[$\times {10}^{-6}{\mathrm{mag}}^{-1}{\mathrm{Mpc}}^{-3}$]
$z\sim 8$
$-22.55$	${0.76}_{-0.41}^{+0.74}$
$-22.05$	${1.38}_{-0.66}^{+1.09}$
$-21.55$aaThis luminosity bin includes sources from the deblending of UVISTA-Y3, which fall below our nominal detection threshold. The sample in this luminosity bin is therefore likely incomplete.	${4.87}_{-1.41}^{+2.01}$
$z\sim 9$
$-22.35$bbThe volume density in this luminosity bin was obtained forcing UVISTA-Y5 into the sample of galaxies at $z_{\mathrm{phot}}\ge 8.6$ (i.e., our $z\sim 9$ LBG sample). Its nominal $z_{\mathrm{phot}}=8.596$ would exclude it from the sample of $z\sim 9$ sources when the redshift selection criteria is strictly enforced; however, considering the very small difference with the $z_{\mathrm{phot}}$ threshold, here we include it for completeness.	${0.43}_{-0.36}^{+0.99}$
$-22.00$	${0.43}_{-0.36}^{+0.98}$
$-21.60$	${1.14}_{-0.73}^{+1.50}$
$-21.20$aaThis luminosity bin includes sources from the deblending of UVISTA-Y3, which fall below our nominal detection threshold. The sample in this luminosity bin is therefore likely incomplete.	${1.64}_{-1.06}^{+2.16}$

Table 7. Table 7: Flux density measurements for UVISTA-Y1 to UVISTA-Y5

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

Table 8. Table 8: Flux density measurements for UVISTA-Y6 to UVISTA-Y12

Filter	UVISTA-Y6		UVISTA-Y7		UVISTA-Y8		UVISTA-Y9		UVISTA-Y10		UVISTA-Y11		UVISTA-Y12
name	[nJy]		[nJy]		[nJy]		[nJy]		[nJy]		[nJy]		[nJy]
CFHTLS $u^{\ast }$	$8$ $\pm$	$14$	$-10$ $\pm$	$15$	$-15$ $\pm$	$15$	$\cdots$		$-10$ $\pm$	$14$	$-3$ $\pm$	$16$	$\cdots$
SSC $B$	$-11$ $\pm$	$9$	$-12$ $\pm$	$9$	$0$ $\pm$	$8$	$4$ $\pm$	$10$	$0$ $\pm$	$9$	$3$ $\pm$	$7$	$-12$ $\pm$	$9$
HSC $g$	$4$ $\pm$	$16$	$-13$ $\pm$	$19$	$0$ $\pm$	$17$	$-2$ $\pm$	$17$	$-3$ $\pm$	$17$	$-19$ $\pm$	$15$	$1$ $\pm$	$16$
CFHTLS $g$	$1$ $\pm$	$13$	$-17$ $\pm$	$13$	$-5$ $\pm$	$13$	$\cdots$		$-21$ $\pm$	$13$	$6$ $\pm$	$12$	$\cdots$
SSC $V$	$0$ $\pm$	$21$	$11$ $\pm$	$21$	$2$ $\pm$	$19$	$-1$ $\pm$	$20$	$11$ $\pm$	$19$	$-1$ $\pm$	$17$	$-1$ $\pm$	$20$
HSC $r$	$11$ $\pm$	$15$	$2$ $\pm$	$21$	$16$ $\pm$	$17$	$3$ $\pm$	$16$	$9$ $\pm$	$16$	$-3$ $\pm$	$14$	$-13$ $\pm$	$16$
CFHTLS $r$	$8$ $\pm$	$21$	$-13$ $\pm$	$23$	$-5$ $\pm$	$19$	$\cdots$		$5$ $\pm$	$19$	$1$ $\pm$	$18$	$\cdots$
SSC $r^{+}$	$20$ $\pm$	$19$	$-34$ $\pm$	$20$	$-21$ $\pm$	$18$	$-8$ $\pm$	$19$	$-14$ $\pm$	$18$	$-1$ $\pm$	$17$	$-18$ $\pm$	$20$
SSC $i^{+}$	$2$ $\pm$	$21$	$-42$ $\pm$	$32$	$-26$ $\pm$	$29$	$-24$ $\pm$	$30$	$-31$ $\pm$	$29$	$-29$ $\pm$	$27$	$-2$ $\pm$	$31$
CFHTLS $y$	$22$ $\pm$	$26$	$-14$ $\pm$	$25$	$-5$ $\pm$	$25$	$\cdots$		$-5$ $\pm$	$25$	$14$ $\pm$	$22$	$\cdots$
CFHTLS $i$	$5$ $\pm$	$29$	$-16$ $\pm$	$27$	$-11$ $\pm$	$26$	$\cdots$		$12$ $\pm$	$25$	$8$ $\pm$	$22$	$\cdots$
HSC $i$	$1$ $\pm$	$23$	$8$ $\pm$	$30$	$-7$ $\pm$	$24$	$-3$ $\pm$	$24$	$-18$ $\pm$	$24$	$14$ $\pm$	$21$	$-16$ $\pm$	$24$
CFHTLS $z$	$-13$ $\pm$	$59$	$-18$ $\pm$	$55$	$-5$ $\pm$	$55$	$\cdots$		$27$ $\pm$	$54$	$-13$ $\pm$	$46$	$\cdots$
HSC $z$	$17$ $\pm$	$33$	$10$ $\pm$	$40$	$-16$ $\pm$	$35$	$-21$ $\pm$	$35$	$50$ $\pm$	$34$	$43$ $\pm$	$32$	$14$ $\pm$	$34$
SSC $z^{+}$	$69$ $\pm$	$85$	$19$ $\pm$	$77$	$16$ $\pm$	$73$	$-74$ $\pm$	$78$	$54$ $\pm$	$72$	$-75$ $\pm$	$63$	$-41$ $\pm$	$76$
HSC $y$	$89$ $\pm$	$80$	$62$ $\pm$	$88$	$55$ $\pm$	$85$	$33$ $\pm$	$87$	$-8$ $\pm$	$85$	$91$ $\pm$	$83$	$14$ $\pm$	$84$
UVISTA $Y$	$16$ $\pm$	$51$	$-9$ $\pm$	$57$	$-31$ $\pm$	$54$	$101$ $\pm$	$53$	$20$ $\pm$	$53$	$-17$ $\pm$	$45$	$-79$ $\pm$	$53$
UVISTA $J$	$211$ $\pm$	$53$	$189$ $\pm$	$64$	$192$ $\pm$	$61$	$204$ $\pm$	$59$	$242$ $\pm$	$61$	$158$ $\pm$	$55$	$164$ $\pm$	$58$
UVISTA $H$	$280$ $\pm$	$66$	$231$ $\pm$	$76$	$241$ $\pm$	$77$	$246$ $\pm$	$78$	$275$ $\pm$	$79$	$297$ $\pm$	$71$	$207$ $\pm$	$76$
UVISTA $K_S$	$271$ $\pm$	$82$	$204$ $\pm$	$84$	$85$ $\pm$	$87$	$116$ $\pm$	$87$	$211$ $\pm$	$88$	$201$ $\pm$	$77$	$191$ $\pm$	$83$
IRAC $3.6\mu$m	$434$ $\pm$	$106$	$235$ $\pm$	$107$	$118$ $\pm$	$81$	$205$ $\pm$	$79$	$153$ $\pm$	$89$	$180$ $\pm$	$54$	$160$ $\pm$	$86$
IRAC $4.5\mu$m	$598$ $\pm$	$130$	$275$ $\pm$	$102$	$310$ $\pm$	$80$	$92$ $\pm$	$100$	$363$ $\pm$	$78$	$376$ $\pm$	$67$	$186$ $\pm$	$88$
IRAC $5.8\mu$m	$-643$ $\pm$	$3000$	$488$ $\pm$	$3604$	$-1188$ $\pm$	$2445$	$1378$ $\pm$	$2596$	$368$ $\pm$	$2423$	$-699$ $\pm$	$3257$	$1815$ $\pm$	$2016$
IRAC $8.0\mu$m	$-3325$ $\pm$	$3803$	$3218$ $\pm$	$5706$	$2049$ $\pm$	$4168$	$-4049$ $\pm$	$2843$	$-207$ $\pm$	$3771$	$2933$ $\pm$	$3972$	$1665$ $\pm$	$3490$

Table 9. Table 9: Flux density measurements for UVISTA-Y13 to UVISTA-Y16

Filter	UVISTA-Y13		UVISTA-Y14		UVISTA-Y15		UVISTA-Y16
name	[nJy]		[nJy]		[nJy]		[nJy]
CFHTLS $u^{\ast }$	$-6$ $\pm$	$13$	$-23$ $\pm$	$18$	$\cdots$		$-28$ $\pm$	$18$
SSC $B$	$-2$ $\pm$	$6$	$-4$ $\pm$	$7$	$8$ $\pm$	$9$	$-6$ $\pm$	$11$
HSC $g$	$-10$ $\pm$	$14$	$4$ $\pm$	$17$	$-7$ $\pm$	$19$	$13$ $\pm$	$16$
CFHTLS $g$	$-2$ $\pm$	$10$	$-1$ $\pm$	$15$	$\cdots$		$-4$ $\pm$	$15$
SSC $V$	$-11$ $\pm$	$17$	$11$ $\pm$	$20$	$-5$ $\pm$	$21$	$18$ $\pm$	$24$
HSC $r$	$5$ $\pm$	$11$	$12$ $\pm$	$13$	$-18$ $\pm$	$17$	$4$ $\pm$	$16$
CFHTLS $r$	$4$ $\pm$	$15$	$2$ $\pm$	$19$	$\cdots$		$-4$ $\pm$	$22$
SSC $r^{+}$	$3$ $\pm$	$15$	$4$ $\pm$	$18$	$-11$ $\pm$	$19$	$-18$ $\pm$	$23$
SSC $i^{+}$	$3$ $\pm$	$24$	$-1$ $\pm$	$28$	$22$ $\pm$	$29$	$10$ $\pm$	$35$
CFHTLS $y$	$14$ $\pm$	$18$	$23$ $\pm$	$26$	$\cdots$		$-16$ $\pm$	$27$
CFHTLS $i$	$-1$ $\pm$	$17$	$8$ $\pm$	$24$	$\cdots$		$1$ $\pm$	$31$
HSC $i$	$0$ $\pm$	$18$	$20$ $\pm$	$22$	$26$ $\pm$	$26$	$-7$ $\pm$	$24$
CFHTLS $z$	$-17$ $\pm$	$43$	$54$ $\pm$	$58$	$\cdots$		$-11$ $\pm$	$64$
HSC $z$	$14$ $\pm$	$26$	$17$ $\pm$	$27$	$47$ $\pm$	$38$	$43$ $\pm$	$35$
SSC $z^{+}$	$-21$ $\pm$	$60$	$29$ $\pm$	$68$	$-41$ $\pm$	$81$	$14$ $\pm$	$75$
HSC $y$	$-4$ $\pm$	$69$	$21$ $\pm$	$81$	$24$ $\pm$	$90$	$-25$ $\pm$	$86$
UVISTA $Y$	$38$ $\pm$	$45$	$-37$ $\pm$	$50$	$69$ $\pm$	$56$	$91$ $\pm$	$56$
UVISTA $J$	$123$ $\pm$	$47$	$144$ $\pm$	$62$	$189$ $\pm$	$56$	$281$ $\pm$	$64$
UVISTA $H$	$166$ $\pm$	$65$	$206$ $\pm$	$75$	$215$ $\pm$	$73$	$268$ $\pm$	$81$
UVISTA $K_S$	$245$ $\pm$	$66$	$34$ $\pm$	$77$	$106$ $\pm$	$77$	$259$ $\pm$	$93$
IRAC $3.6\mu$m	$177$ $\pm$	$96$	$191$ $\pm$	$69$	$134$ $\pm$	$70$	$221$ $\pm$	$79$
IRAC $4.5\mu$m	$362$ $\pm$	$94$	$210$ $\pm$	$93$	$-7$ $\pm$	$83$	$389$ $\pm$	$82$
IRAC $5.8\mu$m	$-287$ $\pm$	$3140$	$-1202$ $\pm$	$3568$	$1349$ $\pm$	$2446$	$-2535$ $\pm$	$2598$
IRAC $8.0\mu$m	$7391$ $\pm$	$4246$	$1553$ $\pm$	$4921$	$-1814$ $\pm$	$3403$	$2196$ $\pm$	$3371$

Table 10. Table 10: V max determination of the UV LF when UVISTA-Y3 is treated as a single source.

$M_{\mathrm{UV}}$	$\varphi$
[mag]	[$\times {10}^{-3}{\mathrm{mag}}^{-1}{\mathrm{Mpc}}^{-3}$]
$-22.55$	${0.0008}_{-0.0004}^{+0.0007}$
$-22.05$	${0.0017}_{-0.0007}^{+0.0012}$
$-21.55$	${0.0036}_{-0.0012}^{+0.0018}$

Table 11. Table 11: V max determination of the UV LF when accounting for lensing magnification.

$M_{\mathrm{UV}}$	$\varphi$
[mag]	[$\times {10}^{-3}{\mathrm{mag}}^{-1}{\mathrm{Mpc}}^{-3}$]
$-22.55$	${0.0008}_{-0.0004}^{+0.0007}$
$-22.05$	${0.0010}_{-0.0006}^{+0.0010}$
$-21.55$	${0.0037}_{-0.0013}^{+0.0019}$

Equations

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

$$J-H>0.8$$

$$(H-K<0.7)\wedge ((K-[3.6]<1.75)\vee (H-[3.6]<1.75))$$

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

Full text

galaxies: formation, galaxies: evolution, galaxies: high-redshift

1 Introduction

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

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

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

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

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

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

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

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

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

2 Observational Data

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

The optical data consists of CFHT/Megacam in $g$, $r$, $i$, $y$ and $z$ (Erben et al. 2009; Hildebrandt et al. 2009 from the Canada-France-Hawaii Legacy Survey (CFHTLS), Subaru/Suprime-Cam $B_{j}$,$V_{j}$, $g^{+}$, $r^{+}$, $i^{+}$ and $z^{+}$-imaging  (Taniguchi et al., 2007), and Subaru HyperSuprimeCam $g$, $r$, $i$, $z$ and $y$ (Aihara et al. 2017a, b).

For this work, we used full-depth Spitzer/IRAC $3.6\mu$m and $4.5\mu$m mosaics we built combining observations from all available programs: S-COSMOS (Sanders et al. 2007), the Spitzer Extended Deep Survey (Ashby et al. 2013), the Spitzer-Cosmic Assembly Near-Infrared Deep Extragalactic Survey (S-CANDELS, Ashby et al. 2015), the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH, PI: Capak), the Spitzer Matching survey of the UltraVISTA ultra-deep Stripes (SMUVS, Caputi et al. 2017; Ashby et al. 2018). Compared to the original S-COSMOS IRAC data, SPLASH provides a large improvement in depth over nearly the whole UltraVISTA area, covering the central 1.2 square degree COSMOS field to 25.5 mag (AB) at 3.6 and $4.5\micron$. SEDS and S-CANDELS cover smaller areas to even deeper limits, while SMUVS pushes deeper over the ultradeep UltraVISTA stripes.

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

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

3 Photometry

Source catalogs were constructed using SExtractor v2.19.5 (Bertin & Arnouts, 1996), run in dual image mode, with source detection performed on the square root of a $\chi^{2}$ image (Szalay et al. 1999) built from the combination of the UltraVISTA $J$, $H$ and $K_{\mathrm{s}}$ images.

The first selection was performed adopting ground-based observations only. Images were first convolved to the $J$-band point-spread function and carefully registered against the detection image (mean RMS $\sim 0\farcs 05$). Initial color measurements were made in small Kron (1980)-like apertures (SExtractor AUTO and Kron factor 1.2) with typical radius $r_{\mathrm{color}}\sim 0\farcs 35-0\farcs 50$.

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

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

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

4 Sample selection

We require sources to be detected at $>5\sigma$ significance in the $J$, $H$, $K_{\mathrm{s}}$, $[3.6]$, and $[4.5]$ images after coadding their S/N’s in quadrature and in those bands with a positive flux density estimate, and we limit our selection to sources brighter than $H\sim 25.8$ mag. The combined UltraVISTA and IRAC detection and S/N requirements exclude spurious sources due to noise, detector artifacts, and diffraction features.

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

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

[TABLE]

while for the $z\sim 9$ sample we required that:

[TABLE]

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

These cuts do not exclusively select $z>7$ galaxies, but also accept some dust-reddened low redshift galaxies. However, such sources would show a very red continuum and red colors red-ward of the $J-$band or $H-$bands. Therefore, to reject this class of galaxies we also imposed to each one of the sample selected with Equations 1 and 2 the requirement of a blue continuuum redward of the break:

[TABLE]

where $\land$ denotes the logical AND operator, and $\lor$ denotes the logical OR operator. These limits are valuable for excluding a small number of very red sources from our selection. Nevertheless, it is worth emphasizing that our final sample of $z>7$ bright galaxies shows little dependence on the specific limits chosen here. This initial selection resulted in 2234 candidates (out of $\sim 8\times 10^{5}$ detected sources): 2015 $Y$ dropouts, and 183 $J$ dropouts.

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

Subsequently, we determined the redshift probability distribution $P(z)$. For this we used the EAzY program (Brammer et al. 2008), which fits non-negative linear combination of galaxy spectral templates to the observed spectral energy distribution (SED), assuming a flat prior on redshifts. We complemented the standard EAzY_v1.0 template set with templates extracted from the Binary Population and Spectral Synthesis code (BPASS - Eldridge et al. 2017) v1.1 for sub-solar metallicity ($Z=0.2Z_{\odot}$), which include nebular emission from cloudy. Specifically, we adopted templates with equivalent widths EW(H$\alpha$)$\sim 1000-5000$ Å as these extreme EW reproduce the observed $[3.6]-[4.5]$ colors for many spectroscopically confirmed $z\sim 7-9$ galaxies (Ono et al. 2012; Finkelstein et al. 2013; Oesch et al. 2015b; Roberts-Borsani et al. 2016; Zitrin et al. 2015; Stark 2016). Driven by current observational results (e.g., Roberts-Borsani et al. 2016; Oesch et al. 2015b; Zitrin et al. 2015), we blanketed the Ly$\alpha$ line from those templates with EW(Ly$\alpha$)$\gtrsim 40$Å. Finally, we added templates of 2 Gyr-old, passively evolving systems from Bruzual & Charlot (2003), with Calzetti et al. (2000) extinction in the range $A_{\mathrm{V}}=0-8$  mag to test the robustness of our selected candidates against being lower-redshift interlopers highly attenuated by dust. We imposed an additional constraint, that the integrated probability beyond $z=6$ to be $>50\%$. The use of a redshift likelihood distribution $P(z)$ is very effective in rejecting faint low-redshift galaxies with a strong Balmer/4000Å break and fairly blue colors redward of the break. After this step, the sample resulted composed of 49 candidates (44 $Y$ dropouts and 5 $J$ dropouts).

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

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

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

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

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

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

5 Results

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

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

5.1 High-resolution imaging from HST

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

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

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

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

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

Given that there are 16 $z\sim 8$ candidates over the $\sim$0.8 deg2 of the UltraVISTA ultradeep stripes, we would expect to find only $\sim$1 candidate over the $\sim$190 arcmin2 CANDELS COSMOS field. Indeed, only one $z\sim 8$ candidate from our selection is located over the CANDELS COSMOS field (UVISTA-Y11). In Figure 5 we present the image stamps in the $V_{606}$ $I_{814}$, $J_{125}$, $JH_{140}$ and $H_{160}$. The $V_{606}$ mosaic shows a close low-z neighbour just $\sim 0\farcs 7$ west of UVISTA-Y11, which is not detected in any NIR image (see Figure 5 and Figure 6). Therefore, we manually included this low-z neighbour when performing the photometry111Omitting the neighbouring source leads to flux densities systematically over-estimated by $\sim 30$%.. We do not detect flux at $>1\sigma$ in the $V_{606}$ and $I_{814}$ bands increasing our confidence on its high-z nature.

Finally, we inspected the ACS $I_{814}$-band mosaic of the COSMOS program (Scoville et al. 2007, $\sim 26.5$ mag in $0\farcs 6$ aperture diameter, $5\sigma$). We found coverage for all sources with the exception of UVISTA-Y1 and UVISTA-Y15. No significant detections exist for any of the sources. We identified a potential low-z galaxy $\sim 1\farcs 0$ north-west of the nominal location of UVISTA-Y4, which however does not affect our flux density estimates.

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

5.2 Sample of $z\sim 8$ Candidates

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

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

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

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

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

5.3 Sample of $z\sim 9$ Candidates

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

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

5.4 Expected Contamination in our Bright $z\sim 8$ Samples

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

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

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

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

5.5 Possible Lensing Magnification

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

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

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

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

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

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

UVISTA-Y6: This source is estimated to be amplified by $\sim$1.4$\times$, $\sim$1.16$\times$ and $\sim$1.14$\times$ from a $10^{10.7}$ $M_{\odot}$, $z=1.76$ galaxy (10:00:12.51, 02:02:57.3), $10^{10.6}$ $M_{\odot}$, $z=1.6$ galaxy (10:00:12.15, 02:02:59.6) and a $10^{10.3}$ $M_{\odot}$, $z=1.65$ galaxy (10:00:12.18, 02:03:00.7), respectively, that lie within $4\farcs 9$, $3\farcs 2$ and $5\farcs 4$ of this source. Their velocity dispersions are estimated to be 259 km/s, 225 km/s, and 206 km/s, respectively.

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

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

UVISTA-Y13: This source is estimated to be amplified by 1.6$\times$ by a $10^{11.15}$ $M_{\odot}$ (330 km/s), $z=1.63$ galaxy (09:58:45.83,01:53:40.6) that lies within $4\farcs 2$ of the source.

We discuss the potential impact of lensing on our inferred value for the characteristic magnitude of the UV luminosity function, $M^{*}$, at the end of Sect. 6.7.

6 Discussion

6.1 Bright candidate LBGs at $z\sim 8$

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

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

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

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

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

Recently, Oesch et al. (2013) analyzed the rest-frame UV and optical properties of a sample of $z\sim 4$ LBGs selected from the GOODS-N/S and HUDF fields and spanning a wide range of UV luminosities, $M_{\mathrm{UV}}\sim-18$ to $\sim-22$ AB. Their $J_{125}-[4.5]$ color (corresponding to approximately rest-frame $u-z$ at $z\sim 4$) shows a correlation with the UV slope $\beta$ (see e.g., their Figure 4), likely driven by dust extinction. The uniform scatter observed at $z\sim 8$ then may suggest rapidly evolving physical mechanisms responsible for the production of dust during the $\sim 800$ Myr between the two epochs.

6.3 Constraints on the EWs of the [OIII]+H$\beta$ lines

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

At redshift $z=7.0-9.1$, the [O III]+H$\beta$ line emission contributes to the Spitzer/IRAC 4.5$\mu$m band in galaxies, producing red $[3.6]-[4.5]$ colors. Figure 13 shows examples of model colors as a function of redshift for lines with very high equivalent width. Using a small sample of $z\sim 8$ galaxies selected from the CANDELS survey, Roberts-Borsani et al. (2016) reported a very red median $[3.6]-[4.5]\sim 0.8$ mag color at bright $H<26$ magnitudes. Using a simple spectral model, consisting of a flat rest-frame $0.3-0.6\mu$m continuum in fν (i.e., a continuum $[3.6]-[4.5]=0$ mag or f${}_{\lambda}\propto\lambda^{-2}$), with the strongest emission lines ([O II]3727, H$\beta$, [O III]4959,5007, H$\alpha$, [N II]6548,6583, [S II]6716,6730), empirical emission lines ratios from Anders & Fritze-v. Alvensleben (2003) for 0.2 Z*⊙* metallicity, they inferred a median [O III]+H$\beta$ EW of $\sim 2000$ Å. However, the sample of Roberts-Borsani et al. (2016) was very small, and possibly biased as it was compiled from IRAC-selected $[3.6]-[4.5]>0.5$ galaxies and galaxies with confirmed Ly$\alpha$ emission. So it is unclear if those results were representative of the general bright $z\sim 8$ population.

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

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

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

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

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

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

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

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

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

We find a positive although marginal correlation of the $M_{\star}/L_{\mathrm{UV,1600}}$ with the $UV$ slope for our $z\sim 8$ sample as it could be expected from older and/or dustier stellar populations characterized by redder UV slope (e.g., Bouwens et al. 2014).

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

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

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

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

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

The total selection volume over our UltraVISTA area for galaxies with $H\sim 24.0-24.5$  mag and $24.5-25.0$ is $5.3\times 10^{6}$ Mpc3 and $2.6\times 10^{6}$ Mpc3, respectively.

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

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

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

For $M_{\mathrm{UV}}\lesssim-22$  mag sources, our new results are also consistent with the upper limits of Bradley et al. (2012), Bouwens et al. (2015), Finkelstein et al. (2015) and of Roberts-Borsani et al. (2016); our measurements are in excess of what is expected extrapolating the Bouwens et al. (2015) results to brighter magnitudes by a factor of $8$, but are nevertheless consistent within $2\sigma$.

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

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

At $M_{\mathrm{UV}}\lesssim-22$  mag our new bright $z\sim 9$ results are consistent with the upper limits of Bouwens et al. (2016) from CANDELS and of Ishigaki et al. (2018) from the Hubble Frontier Field initiative (Lotz et al. 2017). Our measurement at $M_{\mathrm{UV}}\sim-22$  mag is consistent at $1\sigma$ with the measurement of Morishita et al. (2018) based on BoRG observations partly supported by Spitzer/IRAC observations, while it is consistent with that of Calvi et al. (2016) at $\sim 2\sigma$, our density being $\sim 9\times$ lower than the corresponding measurement of Calvi et al. (2016). One possible explanation for differences between our results and those of Calvi et al. (2016) would be if the Calvi et al. (2016) samples suffer from significant contamination. This is especially a concern since few of candidate $z\sim 9$ sources have available Spitzer/IRAC or deep $Y_{098}$ imaging to aid in source selection. In fact, Livermore et al. (2018) find that one especially bright candidate reported by Calvi et al. (2016) appeared to be clearly a low-redshift candidate after further examination.

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

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

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

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

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

In Figure 18 we also compare our estimated Schechter parameters to their evolution over a wide range of redshift, $4\lesssim z\lesssim 10$ from Bouwens et al. (2015). Our result confirm the picture of marginal evolution of $M^{*}$ for $z\gtrsim 3-4$, but significant evolution of $\alpha$ and $\phi^{*}$. This conclusion was first drawn by Bouwens et al. (2015) using LF results from $z\sim 7$ to $z\sim 4$ (see also Finkelstein et al. 2015), although they are in modest tension with the results of Bowler et al. (2015) who suggest an evolution of $dM^{*}/dz\sim 0.2$ from $z\sim 7$ to $z\sim 5$.

6.7 The shape of the LF at $z\sim 8$

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

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

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

[TABLE]

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

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

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

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

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

An alternative way of making sense of the overall shape of the $UV$ LF is to compare it to the halo mass function. To this aim, we scaled the halo mass function by a fixed $M_{\mathrm{halo}}/L_{\mathrm{UV}}$ ratio to match the UV LF. We present the result in Figure 17, adopting the Sheth et al. (2001) halo mass function generated by HMFcalc333http://hmf.icrar.org/hmf_finder/form/create/ (Murray et al., 2013), assuming a Planck Collaboration et al. (2016) cosmology, with $\Omega_{m}=0.2678$, $\Omega_{b}=0.049$, $H_{0}=67.04$, $n_{s}=0.962$, and $\sigma_{8}=0.8347$.

The scaled halo-mass function looks similar to the $UV$ LF at $z\sim 8$. We observe only a small $\Delta\alpha\sim 0.4$ difference in the faint-end slope. We also observe a slight difference in the effective slope at the bright end when a Schechter form is considered ($\Delta\beta\sim 0.2-0.3$). Interestingly, the bright end of the double power-law overlaps with the HMF for $M_{\mathrm{UV}}\lesssim-22$ mag, and might suggest different feedback efficiencies (see also Bowler et al. 2014; Ono et al. 2018). However, as we concluded earlier in this section, our results do not allow us to ascertain whether the cut-off we observe is exponential in form or has a more power-law-like form.

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

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

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

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

In Figure 19 we also present the luminosity density determinations and uncertainties that we derive integrating to $M_{\mathrm{UV}}=-21.5$ mag the recent LF results of Bouwens et al. (2015) at $z\sim 4$, $z\sim 5$, $z\sim 6$, $z\sim 7$ and $z\sim 8$, the results of Bouwens et al. (2016) at $z\sim 9$, and at $z\sim 10$ from Oesch et al. (2018), together with our best-fit constraints and $1\sigma$ uncertainties on the evolution of the luminosity density with redshift, assuming that the logarithm of the luminosity density decreases linearly with increasing the redshift.

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

7 Conclusions

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

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

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

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

We note that, despite much recent progress, the lack of spectra and deep high-resolution imaging still limit us in establishing the reliability of high redshift galaxy selections, in particular for rare luminous galaxies that constrain the bright end of the UV luminosity and mass functions where any contamination has a very large impact. While care is taken in estimating the completeness and contamination rates, these still rely on assumed spectral energy distributions. Ultimately, spectroscopy is needed to validate these assumptions. While recent results suggest ALMA as a potentially efficient machine for the study of emission lines (e.g., Smit et al. 2018), currently this is still hard in the rest-frame UV and optical, due to long integration times, low multiplexing, and the reduced observable emission of Ly$\alpha$ likely caused by the increasing neutral hydrogen fraction $z>6$ (e.g., Schenker et al. 2014), but will be possible in the future with JWST and next generation of extremely large ground-based telescopes.

The authors would like to thank the referee for their careful reading and for the many helpful, positive and very constructive comments that helped improving the quality of the paper. The authors wish to acknowledge the UltraVista team who conducted and followed the observations and image processing. The authors are thankful to the COSMOS collaboration for their continued efforts in making available part of the data used in this work. KIC acknowledges funding from the European Research Council through the award of the Consolidator Grant ID 681627-BUILDUP. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. The Cosmic Dawn Center is funded by the Danish National Research Foundation. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Based in part on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan.The 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.

Appendix A Monte Carlo assessment of the multi-component nature of UVISTA-Y3

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

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

In Figure 20 we present the twenty random realizations of the exponential disk, before convolution with the WFC3 PSF, while in Figure 21 we present the image stamps of the DASH mosaics after the synthetic exponential disks have been added. Based on simple visual inspection, we see no indication in these simulated images for a multiple component structure. Finally we run SExtractor using the same set of parameters adopted for the original deblending and found that none of the synthetic sources were split into two or more components. This test therefore increased our confidence on the multiple nature of UVISTA-Y3. It is worth remarking that as a result of the low S/N of the deblended photometry for each component there are substantially larger uncertainties in the derived physical parameters for each component.

We complemented this first assessment with a second Monte Carlo simulation in which we adopt an effective radius $r_{\mathrm{e}}=3$ kpc, similar to the sizes of luminous high redshift LBGs when potentially multiple sources are considered as a single object (e.g., Bowler et al. 2017), and $0.1<b/a<0.7$. We present the result of this simulation in Figures 22 and 23. Using a procedure similar to that applied for the $r_{\mathrm{e}}=1$ kpc case, we do not find evidence for multiple components even when large $r_{\mathrm{e}}$ are considered, increasing the confidence on our interpretation of the three components in UVISTA-Y3.

Appendix B Flux density estimates

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

Appendix C LF estimate when UVISTA-Y3 is considered as one single source

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

Appendix D LF estimate after correcting for lensing magnification

Here we present the $z\sim 8$ LF determination after correcting for lensing magnification the luminosities of four objects (Y6, Y8, Y9 and Y13 - see Section 5.5). None of the sources affected by lensing satisfies the criteria for inclusion in the $z\sim 9$ LF.

The change of solid angle introduced by the lenses affects the volume estimates and therefore the volume densities. A proper estimate of the volumes that takes into account the lensing effects would require running a simulation with the actual lensing effects from all the sources in the region considered for our search. While this goes beyond the scope of this paper, we can attempt to estimate the average effect based on the statistics in our sample. Considering that $4/18$ objects in our sample have cumulative lensing magnification factors of $\sim 1.85,1.96,1.39$ and $1.60$ respectively, conserving the surface brightness would imply shrinking by the same amounts the area corresponding to each source (i.e., $\sim 1/18$ the area of the UltraVISTA field). Given that the remaining $14/18$ sources are not significantly magnified, the total area available for the search would then be $91\%$ of the area when no lensing is considered. The volumes densities would then be $\sim 10\%$ higher (just $\sim 0.04$ dex) than those computed without introducing the lensing magnification. Because of the very marginal effect of this correction compared to the Poissonian and cosmic variance uncertainties, the LF measurements presented in this Section were computed with the same volumes adopted for the LFs presented in Section 6.5, which do not include any correction for lensing effects. In Table 11 we list the corresponding $V_{\mathrm{max}}$ determination when the luminosities of those sources are corrected for lensing magnification.

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