The Bright-end Galaxy Candidates at z ∼ 9 from 79 Independent HST Fields

BoRG[z9] is a pure-parallel HST imaging survey with five broadband filters spanning from the near-UV/optical band of WFC3/UVIS (${O}_{350}$) to near-IR (${Y}_{105}$/${J}_{125}$/${{JH}}_{140}$/${H}_{160}$) of WFC3/IR. Its medium-deep exposures (∼2–5 ks) in random sightlines optimize for detecting bright galaxies, and reach typical limiting magnitudes of ${m}_{5\sigma }\sim 27$ mag measured in an aperture of 032 radius (see Appendix A).

The major update in BoRG[z9] from our previous studies (HST Cycles 17 and 19; Trenti et al. 2011, 2012a; Bradley et al. 2012; Schmidt et al. 2014) is the use of a long-pass filter in the optical wavelength range. ${O}_{350}$ covers ∼0.3 μm to 1.0 μm, while many other studies use a single or multiple broadband filters in the optical range. The choice optimizes rejection of contaminants for the selection of $z\gt 8$ sources.

For our primary science goal, the collection of ${J}_{125}{H}_{160}$ and ${Y}_{105}{{JH}}_{140}$ dropout sources, persistence is of particular concern because of the possibility of introducing an artificial coherent signal into the near-IR bands. In each visit of this program, we therefore arranged the sequence of WFC3/IR filters to minimize the impact. As detector persistence decays over time, with approximate power-law behavior, any saturated target observed in a previous visit most affects the initial part of the pure-parallel orbit. The general strategy therefore is to observe in the dropout filter as early as possible in the orbit, ideally placing the ${Y}_{105}$ band first. However, doing so would result in increased time-varying backgrounds from scattered-light Earth glow from helium in the upper atmosphere (caused by a He 10830 Å line; Brammer et al. 2014). Because the intensity of the scattered Earth glow decreases with increasing target angles above the bright Earth limb, we opted to place the ${Y}_{105}$ observations as the second filter after ${J}_{125}$ (or sometimes third when the orbit begins with ${O}_{350}$) in an orbit. Whenever possible, observations in ${{JH}}_{140}$ and ${H}_{160}$ follow those filters. With this strategy, persistence features are essentially prevented from contaminating the ${J}_{125}{H}_{160}$ and ${Y}_{105}{{JH}}_{140}$ dropout selections.

To ensure good sampling of the IR exposures, we opted for reading every 50 s (SPARS50). While the majority of cosmic rays are rejected by the calibration pipeline, owing to the multiple non-destructive readouts of the WFC3IR detector, a small fraction may survive in the calibrated image. We thus split the total integration in each IR filter into at least two individual exposures. ${O}_{350}$ exposures obtained with the UVIS CCD are split into more than two sequences (each 450–900 s) for optimal cosmic-ray rejection. Our design choices are aimed at maximizing the data quality, though a small price is paid in signal-to-noise ratio (S/N) by increasing the number of individual exposures.

We reduced the raw imaging data by using the HST pipeline in the standard manner. In addition to running calwf3, processing of ${O}_{350}$ included a correction for the charge-transfer efficiency effect (CTE19; Noeske et al. 2012; Anderson 2014). For all filters, we also performed a customized extra step to remove residual cosmic rays and/or detector artifacts, such as unflagged hot pixels, by using a Laplacian edge filtering algorithm (LACOSMIC; van Dokkum et al. 2001) that was used for the previous BoRG data reduction (Bradley et al. 2012; Schmidt et al. 2014, see also C16).

The pixel scale is set to 0.08 arcsec/pixel with pixfrac to 0.75, as in C16. Since most of our images are taken without dithering, we use a slightly larger pixel scale than typical non-parallel observations (∼0.06 arcsec/pixel). The only change from C16 is made in one of the parameters of cosmic-ray detection, where we set cr_threshold to 3σ and 1.5σ (compare 3.9σ and 5σ in C16) for UVIS and IR detectors, respectively, in the LACOSMIC Python package.¹⁴ This significantly reduces the residual cosmic rays compared to the previous data products.

This reduction process generates science and rms maps for the five filter bands, in addition to a combined ${{JH}}_{140}$ + ${H}_{160}$ map for the use of detection images (see below). It is noted that the rms maps generated by the pipeline have arbitrary infinity values for, e.g., dead pixels and/or those with a cosmic-ray flag. Such an artificial value in fact affected the previous sample selection in C16 (see Section 3.1.1). We replace those values with 0, so that a detection algorithm (i.e., SExtractor; Bertin & Arnouts 1996) ignores them; otherwise it would falsely return infinity values for fluxes in such regions.

The rms map of each filter band generated by the pipeline is then scaled so that it represents the true uncertainty, including correlated noise (Casertano et al. 2000). To account for this, we follow the scaling method presented by Trenti et al. (2011). Briefly, we measure the median sky in the empty region in each mosaic image by using SExtractor (Bertin & Arnouts 1996) with the same setup parameters as the source detection. Each rms map is then scaled so that the dispersion in flux of the sky region (f_aper, in an aperture of $r=0\buildrel{\prime\prime}\over{.} 32$) and its estimated error (e_aper) are consistent. The median values of the scale factor are 1.174, 1.195, 1.209, 1.233, and 1.708 for ${H}_{160}$, ${{JH}}_{140}$, ${J}_{125}$, ${Y}_{105}$, and ${O}_{350}$, respectively, which are consistent with C16. This process also returns limiting magnitudes (i.e., median of f_aper in each field), which are used for upper limits of non-detection.

Ten out of the 89 original survey fields are discarded from our final analysis, because they failed in the acquisition of a guide star or are dominated by a number of bright stars. The list of fields is summarized in Appendix A.

The point-spread functions (PSFs) of different filter images are not matched in this study, because it is challenging to accurately measure the PSF in each of our data with relatively shallow exposures. However, as the apparent size of our target galaxies is typically small, in addition to the use of optimized photometry (i.e., isophotal flux; see Section 2.3), the effect is minimal (see also Oesch et al. 2007; Trenti et al. 2012a; Holwerda et al. 2014).

We first detect galaxies using SExtractor in the ${{JH}}_{140}$ + ${H}_{160}$ stacked image. The detection parameters are set as in C16—DETECT_MINAREA = 9, NTHRESH = 0.7σ, DEBLEND_NTHRESH = 32, DEBLEND_MINCONT = 0.01. The exception is the convolution size (FWHM of Gaussian) for detection in SExtractor, which is changed from 2 pixels to 5 pixels. This reduces false detection of, e.g., discrete noise at the edge of the detector and residual cosmic rays. The photometry is performed in the dual-imaging mode, based on the ${{JH}}_{140}$ + ${H}_{160}$ detection. With this setup, we detect 73,374 objects from all BoRG fields in this study.

We then select those with ${\rm{S}}/{{\rm{N}}}_{{{JH}}_{140}+{H}_{160}}\gt 8$ for robust detection, and with CLASS_STAR < 0.95 to avoid stars. Signal-to-noise ratios are calculated from the measurement and error of the isophotal flux of SExtractor. Compared to other schemes of photometry, the choice flexibly corresponds to source morphology. For example, fixed-aperture photometry tends to underestimate S/N in ${O}_{350}$, which would mistakenly select low-z galaxies as the dropout candidates.

To minimize contamination by low-z interlopers, we also limit half-light radius along the major axis of SExtractor to 03 (or ∼1.7 kpc at $z\sim 6$), the observed upper limit for $z\gt 6$ galaxies at the present magnitude limit (e.g., Oesch et al. 2010; Holwerda et al. 2015; Shibuya et al. 2015; Kawamata et al. 2018).

We then follow the color-cut criteria proposed by C16 for the selection of z ∼ 10 and z ∼ 9 candidates, whose 95% confidence regions correspond to $7.7\lt z\lt 9.7$ (with a peak at $z\sim 8.7$) and $z\gt 9.3$, respectively;

• 1.  
  z ∼ 10 candidates (${J}_{125}{H}_{160}$ dropouts):

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{350}\lt 1.5,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{105}\lt 1.5,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{160}\geqslant 6,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{J}_{125}-{H}_{160}\gt 1.3.\end{eqnarray*}$$
• 2.  
  z ∼ 9 candidates (${Y}_{105}{{JH}}_{140}$ dropouts):

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{350}\lt 1.5,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{140}\geqslant 6,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{\rm{S}}/{{\rm{N}}}_{160}\geqslant 4,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{Y}_{105}-{{JH}}_{140}\gt 1.5,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{Y}_{105}-{{JH}}_{140}\gt 5.33({{JH}}_{140}-{H}_{160})+0.7,\end{eqnarray*}$$

  $$\begin{eqnarray*}&&{{JH}}_{140}-{H}_{160}\lt 0.3.\end{eqnarray*}$$

Colors are calculated after correcting for Galactic dust extinction, where E(B − V) is retrieved from the NASA/IPAC infrared archive (Schlegel et al. 1998; Schlafly & Finkbeiner 2011).¹⁵ We assume the canonical Milky Way dust law (Cardelli et al. 1989) to calculate the extinction in each filter.

The observed flux is scaled by ${C}_{160}={f}_{\mathrm{AUTO},160}/{f}_{\mathrm{iso},160}$ for the following analysis, where ${f}_{\mathrm{AUTO},160}$ and ${f}_{\mathrm{iso},160}$ are the AUTO (i.e., total) flux and isophotal flux of SExtractor measured in ${H}_{160}$, so as to correct the flux missed in the isophotal flux. The scale factor of individual objects is applied to the other four bands to uniformly scale the fluxes.

The selected sources are shown in Figure 1. There are 15 and 19 sources selected with the color criteria for z ∼ 10 and z ∼ 9 candidates, respectively, up to ${H}_{160}\sim 27$ mag. The color near the Lyman break (${J}_{125}{H}_{160}$ and ${Y}_{105}{{JH}}_{140}$) is calculated with the 1σ limiting magnitude of each image when the magnitude at the shorter wavelength has ${\rm{S}}/{\rm{N}}\lt 1$, and shown as a lower limit in Figure 1.

Zoom In Zoom Out Reset image size

While our color selection mostly excludes sources with optical detection, we found in the visual inspection (Section 2.4) that some candidates show tiny blobs in ${O}_{350}$, despite satisfying ${\rm{S}}/{{\rm{N}}}_{350}\lt 1.5$ in the isophotal magnitude. Given its bright/extended appearance in longer wavelength filters, we consider this as sub-galactic scale/patchy star formation from dusty galaxies at lower redshifts (see Section 2.5). We therefore exclude those with ${\rm{S}}/{{\rm{N}}}_{350,\mathrm{ap}.}\gt 1.5$ in addition to the criteria above, where ${\rm{S}}/{{\rm{N}}}_{350,\mathrm{ap}.}$ is a signal-to-noise ratio measured in a small aperture (radius of 2 pixels, or 016). The aperture photometry maximizes S/N inside objects for this case, while isophotal flux is typically measured in a larger aperture from detection images (${{JH}}_{140}$ + ${H}_{160}$ band in this study) and could lower ${\rm{S}}/\mathrm{Ns}$ at shorter wavelength. This excludes one of the z ∼ 9 candidates selected with the default criteria, while there are no such objects in z ∼ 10 candidates. The remaining candidates are taken through the following selection process.

Four of the authors (T.M., M.T., M.S., and R.C.L.) visually checked the color-selected candidates, because we are not able to exclude artificial objects, such as PSF spikes and residual cosmic rays with only the color criteria. We exclude five objects selected as z ∼ 10 candidates and five selected as z ∼ 9 candidates. Those are stellar spikes, cosmic rays, or inappropriately deblended fragmentation of large galaxies.

We also check the persistence among candidates by inspecting the individual fields.¹⁶ We conclude that none of the remaining candidates is affected by persistence.

2.5.1. A Flat Prior

2.5.2. Building a Practical Prior

A flat prior is often assumed for Lyman break galaxies at high z, since it is still unclear whether applying priors from low redshift is appropriate (Benítez 2000; Salmon et al. 2017; Salvato et al. 2018). On the other hand, we know that the LF of galaxies evolves strongly from z ∼ 10 to 0 (e.g., Behroozi et al. 2013; Bouwens et al. 2015b; Mason et al. 2015a; Kelson et al. 2016; Oesch et al. 2018). Combining the intrinsic evolution with changes in the distance modulus with redshift, we should, at a given luminosity, expect many more galaxies at lower redshift in a given survey area (Stiavelli 2009; Vulcani et al. 2017). This effect has not been taken into account yet in this study.

Motivated by this fact, we design an empirical prior based on a combination of deep observations and theoretical modeling, which aims at quantifying the relative abundance of high-z sources versus low-z interlopers for color-selected samples using the dropout technique. To construct such a prior, we start from a mock galaxy catalog by Williams et al. (2018). While the catalog (v1.0) has been primarily aimed at planning observations with the James Webb Space Telescope (JWST), based on real deep surveys with multiband imaging including most HST broadband filters, it reproduces observed galaxy properties self-consistently, such as LF evolution up to $z\sim 8$. Therefore, it is well suited for our purposes.

The model flux in the catalog for each galaxy is perturbed to simulate observed fluxes. We assume Gaussian noise based on our limiting magnitudes. In addition to this random noise added to the model flux, we fluctuate the intrinsic flux by 0.3 dex at each redshift bin (whose width is δz = 0.1), to partially take into account the effect of cosmic variance (Somerville et al. 2004; Trenti & Stiavelli 2008), which is not reflected in the original catalog based on relatively small-volume observations. While this is a simplified approach, shifting magnitude artificially adds ∼0.4 dex variation in galaxy number densities at ${M}_{\mathrm{UV}}\lesssim -20$. We repeat this 10,000 times (i.e., 10,000 mock catalogs with different cosmic fluctuations, where each object has Gaussian random noise added to its flux), then select dropout candidates with our color criteria from all catalogs. The only difference is the use of F435/606/814W bands in the mock catalog, because ${O}_{350}$ flux is not listed in the original catalog. We estimate the weighted mean flux from the three optical bands as

$$\begin{eqnarray}&&{f}_{350}=\displaystyle \sum _{i}{f}_{i}\times {w}_{i}/\displaystyle \sum _{i}{w}_{i}\end{eqnarray} \tag{ 1 }$$

where the weight is the convolution of transmission in each filter with ${O}_{350}$,

$$\begin{eqnarray}&&{w}_{i}={T}_{i}\ast {T}_{350}.\end{eqnarray} \tag{ 2 }$$

From the mock catalog, we estimate the fraction of low-z interlopers that are selected with our color selection, as a function of magnitude and redshift. We define low-redshift interlopers as those with intrinsic redshift at $z\lt 8$ but which satisfy the color–color selection criteria in Section 2.3.

Figure 2 shows examples—the intrinsic redshift distributions of the dropout candidates at a given observed ${H}_{160}$ magnitude. The fraction of interlopers ranges from ∼70% (∼40%) to ≳85% (∼90%) for z ∼ 10 (z ∼ 9) candidates from ${H}_{160}$ ∼ 24.8 mag to 26.5 mag, which is consistent with the analysis of Vulcani et al. (2017).

Zoom In Zoom Out Reset image size

The interloper population is dominated by z ∼ 2–3 galaxies. The fraction of star-forming/quiescent populations in the low-z interlopers is shown in the bottom panel of Figure 2. While the low-z interlopers of z ∼ 10 candidates consist of both populations, i.e., the star-forming population at $z\sim 2$ and quiescent population at $z\sim 3.5$, those of z ∼ 9 candidates are dominated by the star-forming population at $z\sim 3$. For ${Y}_{105}{{JH}}_{140}$ dropout candidates, the fraction of low-z interlopers increases as it approaches the magnitude ${H}_{160}$ ∼ 26.5, because large photometric flux errors start to affect the sample selection. On the other hand, for ${J}_{125}{H}_{160}$ dropout candidates it shows a rather flat distribution over the magnitude range. This is partly due to different selection effects for the color selection criteria—while ${Y}_{105}{{JH}}_{140}$ dropout candidates are mainly affected by photometric error, ${J}_{125}{H}_{160}$ dropout candidates are contaminated mainly by systematic effects due to, e.g., a smaller number of filters to characterize the Lyman break.

2.5.3. Application of the Redshift Prior

In general, a redshift posterior is calculated as

$$\begin{eqnarray}&&p(z| C,m)\propto p(z| m)p(C| z),\end{eqnarray} \tag{ 3 }$$

where $p(C| z)$ is a likelihood (derived by EAzY) given the data C. $p(z| m)$ is a prior as a function of m (e.g., magnitude; Benítez 2000). In the present study, since we first select galaxies with color criteria, an additional term is implicitly implemented in Equation (3):

$$\begin{eqnarray*}p(C)=\left\{\begin{array}{ll}1, & \mathrm{if}\ \mathrm{objects}\ \mathrm{satisfy}\ \mathrm{the}\ \mathrm{color}\ \mathrm{selection}\\ 0, & \mathrm{otherwise}.\end{array}\right.\end{eqnarray*}$$

With this, redshift distributions derived in the previous section (Figure 2) can be used as priors, $p(C)p(z| m)$, in the calculation of posterior for our preselected candidates.

To build a functional form for priors, we fit the redshift distribution with a dual-peak Gaussian model, $p(C)p(z| m)={\sum }_{i=\mathrm{1,2}}{G}_{i}(z| {I}_{i},m)$, where ${I}_{i}=\{{A}_{z,i},{\mu }_{z,i},{\sigma }_{z,i}\}$ is the parameter set for the Gaussian fit, ${G}_{i}(z| {I}_{i})={A}_{z,i}\,\exp [-{(z-{\mu }_{z,i})}^{2}/2{\sigma }_{z,i}^{2}]$. The best-fit parameters (${\chi }^{2}$ minimization) are summarized in Table 1. Due to the small number of galaxies, Gaussian fits for z ∼ 10 candidates in the ${H}_{160}\lesssim 25$ and ≳26.3 mag bins are not available. We instead extrapolate without changing the parameter values from the 25 mag and 26.25 mag bins to those magnitude bins, respectively.

Table 1.  Gaussian Parameters for Redshift Priors

m160	${A}_{z,1}$	${\mu }_{z,1}$	${\sigma }_{z,1}$	${A}_{z,2}$	${\mu }_{z,2}$	${\sigma }_{z,2}$	${f}_{\mathrm{int}.}$ a
${J}_{125}{H}_{160}$ dropouts
25.00	0.87	3.29	0.57	0.13	10.43	0.58	0.83
25.25	0.66	3.13	0.60	0.34	10.53	0.12	0.86
25.50	0.57	3.13	0.62	0.43	10.51	0.12	0.82
25.75	0.74	2.90	0.56	0.26	10.24	0.52	0.76
26.00	0.77	2.51	0.61	0.23	10.11	0.62	0.77
26.25	0.89	2.23	0.43	0.11	9.97	0.69	0.86
${Y}_{105}{{JH}}_{140}$ dropouts
24.50	0.43	2.12	0.08	0.57	8.76	0.14	0.30
24.75	0.56	2.09	0.10	0.44	8.72	0.16	0.37
25.00	0.57	2.07	0.10	0.43	8.72	0.17	0.38
25.25	0.61	2.06	0.09	0.39	8.75	0.26	0.36
25.50	0.70	2.16	0.06	0.30	8.76	0.34	0.34
25.75	0.55	2.08	0.12	0.45	8.74	0.46	0.35
26.00	0.60	2.08	0.16	0.40	8.70	0.51	0.44
26.25	0.75	2.08	0.20	0.25	8.61	0.57	0.65
26.50	0.86	2.07	0.28	0.14	8.47	0.66	0.84
26.75	0.88	1.99	0.59	0.12	8.31	0.78	0.92

Notes. For ${J}_{125}{H}_{160}$ dropout candidates with ${m}_{160}\lt 25$ and $\gt 26.25$ mag, reliable Gaussian fitting parameters are not available from the small sample size. We instead adopt those from ${m}_{160}=25$ and 26.25 mag for candidates, respectively.

^aFraction of low-z interlopers.

Download table as:  ASCIITypeset image

Given the smaller error in ${H}_{160}$ magnitudes (<0.2 mag) compared to the magnitude bin of the prior, we simply multiply the fitted Gaussian function by each redshift probability distribution derived with a flat prior in the earlier part of this section. We use the redshift at the peak of the posterior distribution as the final redshift (z_prior). We exclude those candidates with a low-z peak probability of >20% of the high-z peak probability, based on the typical contamination fraction at z ∼ 10 (Pirzkal et al. 2013). Thus, the photo-z selection assures 80% reliability, or equivalently a contamination fraction ${f}_{\mathrm{cont}}=0.2$ (Section 4.2). We reject nine and five of the remaining z ∼ 10 and z ∼ 9 candidates, respectively.

While the mock catalog of Williams et al. (2018) reproduces the observed LFs at $z\lesssim 8$ sufficiently well, we note that galaxies at $z\gt 8$, especially at the bright end, are yet uncertain, and an a priori assumption, such as the Schechter-form LF, may no longer be appropriate. It is also possible that samples might be somewhat affected by currently unaccounted for contamination by a class of rare intermediate-redshift galaxies with observed colors similar to those of $z\gtrsim 9$ sources. Yet, all our final candidates have a low likelihood of being interlopers (less than 20%), regardless of the application of the prior, strengthening confidence in their high-z nature. In fact, the application of the prior suppresses high-redshift probabilities for most of our cases (Table 1), and thus it offers a more conservative candidate selection. Of course, further follow-up of our candidates with deeper imaging, or spectroscopy, would represent an independent assessment of the selection.

For further validation of our candidates, we check the availability of public Spitzer data¹⁷ acquired from multiple Spitzer/IRAC 3.6 μm observations of BoRG fields (PIs S. Bernard, R. Bouwens, B. Holwerda; ${m}_{5\sigma }\sim 25\mbox{--}26.5$ mag). While the data set does not offer complete coverage of the entire BoRG[z9] survey, IRAC photometry is capable of efficiently excluding low-z interlopers in the case of significant excess of flux compared to the ${H}_{160}$ band such as, e.g., quiescent+dusty SEDs at $z\sim 2.5$ (Holwerda et al. 2015; Oesch et al. 2016; Salmon et al. 2018).

The z ∼ 10 candidate has an IRAC 3.6 μm coverage. This candidate is located near a bright galaxy, where the photometry in the IRAC band is challenging with aperture photometry (Figure 3). To extract the flux from the candidate, we use TPHOT (Merlin et al. 2016), which models the IRAC flux based on the light profile obtained from the high-resolution HST images. We use the ${{JH}}_{140}$ + ${H}_{160}$ image as the reference model and a convolution kernel that is constructed from the PSFs in the ${H}_{160}$ and IRAC 3.6 μm bands. This approach yields an extracted model magnitude of ${m}_{3.6}=23.8\pm 0.7$ for the candidate. By including this IRAC flux in the photometric redshift code, the low-z probability peak is ∼16% of the high-z peak, and thus we retain the galaxy as the final, and only, z ∼ 10 candidate in this study (Figure 3).

Zoom In Zoom Out Reset image size

One of the z ∼ 9 candidates, 0956+2848–98, has IRAC 3.6 μm coverage, but shows no detection (${\rm{S}}/{\rm{N}}\lt 1$), and a 1σ lower limit is available from the rms map (>25.5 mag; Figure 4). The photometric redshift of the source including the IRAC upper limit remains unchanged, and thus we retain this in the sample of final z ∼ 9 candidates.

Zoom In Zoom Out Reset image size

One of our z ∼ 9 candidates is by chance overlapping with previous archived observations, both from the Hubble Infrared Pure Parallel Imaging Extragalactic Survey (HIPPIES) pure-parallel program (GO 11702; PI: Yan; ∼10 orbits) and from a previous BoRG campaign (∼1 orbit). Combining all data, and in particular thanks to the very deep GO 11702 exposure, the object shows a significant detection in the deep F098M image, ${m}_{098}=27.1\pm 0.2$. The most likely photometric redshift for this object becomes $z\sim 2.2$, and thus we exclude it from the final candidate sample. Two other z ∼ 9 candidates, which were reported in C16, have also been recently followed up with the same filter (1 orbit; Livermore et al. 2018). The follow-up data revised the redshift to $z\sim 1.8$ for one of the two galaxies (while confirming the high-z nature of the other), and thus this source is not included in our final sample.

This suggests that, even after the stringent selection of this study, low-z interlopers could be selected as final high-z candidates. While our final candidates show low probability peaks at low z (≲16%), additional imaging around the putative Lyman break would sample the SED more accurately and further improve the quality of the high-z photometric candidates in the absence of spectroscopic confirmation. Deep Spitzer data would constrain the rest-frame optical wavelength range, but only if sources are well isolated (see the case for 2140+0241–303). Space-based deep medium-/narrowband observations, as in the BoRG Cycle 25 program (Livermore et al. 2018), are ideal before the advent of next-generation facilities such as JWST (see also Appendix C).

The UV absolute magnitude at rest-frame 1450 Å, M_UV, is then calculated based on the photometric redshift (i.e., the k-correction) and scaled ${H}_{160}$ magnitude (Section 2.3). We assume a log-normal SED with a UV slope β = −2.0 (${f}_{\lambda }\propto {\lambda }^{\beta }$), appropriate for the high-z population (e.g., Fan et al. 2003; Bouwens et al. 2014). We checked that our choice of β from observed values (∼−1.5 to −2.0 for ${M}_{\mathrm{UV}}=-23$ to −20; Bouwens et al. 2014) does not change our final result.

Among the sample, three candidates have neighboring foreground objects in each FoV. The observed light from those high-z candidates is thus possibly affected by the gravitational potential of the foreground objects.

We calculate the magnification by foreground objects in the same manner as Mason et al. (2015b). Briefly, the photometric redshift of low-z objects is derived with EAzY. The probability distributions for both foreground and high-z candidates, and their separation, are then used to estimate the magnification. Single isothermal spheres are assumed for the mass profile of foreground deflectors. The Einstein radius of deflectors is estimated from photometry using a redshift-dependent relation of Faber & Jackson (1976). The resulting magnifications are ${1.5}_{-0.3}^{+0.7}$ (2140+0241−303), ${1.1}_{-0.1}^{+0.1}$ (0751+2917−499), and ${1.7}_{-0.4}^{+0.5}$ (2229−0945−394). The error in magnification is also integrated into the calculation of absolute magnitude.

Part of our data (36 fields; ∼40%) has been studied in C16. While the initial color selections are identical in both studies, the application of our selection processes, in addition to the update in drizzle pipeline parameters and SExtractor parameters, may affect the final sample. In what follows, we compare our candidates with those presented in C16.

3.1.1. z ∼ 10 Candidates

3.1.2. z ∼ 9 Candidates

We calculate the effective volume by following Oesch et al. (2012) and C16 (see also Carrasco et al. 2018, for a public code based on the same approach adopted here). Briefly, we added artificial sources in empty regions of each science frame. Those sources were modeled with realistic intrinsic distributions of UV colors and half-light radii, and they had assigned redshifts and ${H}_{160}$ magnitudes to compute the selection efficiency as a function of those parameters.

The sources were modeled with a mix of n = 1 and n = 4 Sérsic (1968) profiles, and with half-light radii as found in previous studies (Oesch et al. 2010; Grazian et al. 2012; Ono et al. 2013; Holwerda et al. 2015; Kawamata et al. 2015; Curtis-Lake et al. 2016; Bouwens et al. 2017). We assumed that there is no significant evolution at $z\gtrsim 8$ (Wilkins et al. 2016), and modeled the SED with observed UV slope β fixed to the $z\sim 8$ value (Bouwens et al. 2014). We calculated the completeness as a function of apparent ${H}_{160}$ magnitude, C(m), and the source selection function as a function of magnitude and redshift, $S(z,m)$, in each sightline by detecting those artificial sources.

The effective comoving volume is then calculated as

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}(m)={\int }_{0}^{\infty }S(z,m)C(m)\displaystyle \frac{{dV}}{{dz}}{dz},\end{eqnarray} \tag{ 4 }$$

where dV/dz represents differential comoving volume at redshift z. The total volume, the sum of the effective volume over all the effective survey fields, i.e., the region that is not occupied by bright sources, ranges from ∼4 × 10⁴ Mpc³ to ∼9 × 10⁵ Mpc³ for our magnitude range, $-24\lesssim {M}_{\mathrm{UV}}\,\lesssim \,-20$. The effective volume used for each magnitude bin is summarized in Table 3.

Table 3.  Number Density of Dropout Candidates

MUV	$\,\mathrm{log}\phi$	Number	${V}_{\mathrm{eff}}$
(mag)	(Mpc−3 mag−1)		(104 Mpc3)
z ∼ 10
−23.0	$-{6.1}_{-0.8}^{+0.5}$	1	95.68
−22.0	<−5.9	0	76.08
−21.0	<−4.6	0	3.80
z ∼ 9
−23.0	<−5.9	0	73.39
−22.0	$-{5.9}_{-0.8}^{+0.5}$	1	63.87
−21.0	$-{5.4}_{-0.8}^{+0.5}$	1	19.94

Note. Errors of number densities are dominated by Poisson uncertainty.

Download table as:  ASCIITypeset image

With the final candidates and effective volume derived in the previous sections, we estimate the stepwise number density. The density is estimated by dividing the number of candidates at a given UV absolute magnitude by the effective volume. The 1σ confidence level is estimated by assuming a Poisson distribution (Gehrels 1986). For magnitude bins with no candidates, we show 1σ upper limits derived from the Poisson distribution and effective volume.

We also take account of contamination by multiplying by $(1-{f}_{\mathrm{cont}})$, where ${f}_{\mathrm{cont}}=0.2$ is the contamination fraction defined in Section 2.5 (see also Bradley et al. 2012; Schmidt et al. 2014). With this, the uncertainty from this arbitrary contamination fraction is now limited to ≲0.1 dex in the estimation of number density.

Our results are summarized in Table 3 and shown in Figure 5, where we find consistency with previous studies in BoRG (C16) and other fields (Oesch et al. 2013, 2018; Bouwens et al. 2015b). Our upper limits constrain the number density at values a factor of ∼0.3 dex lower than in C16, thanks to the volume from 79 independent fields. It is noted that since three of C16's candidates are rejected as low-z interlopers in this study, one of C16's data points at z ∼ 10 (${M}_{\mathrm{UV}}\sim -22.3$) now becomes an upper limit.

Zoom In Zoom Out Reset image size

At z ∼ 9, we find a good agreement with theoretical expectations at ${M}_{\mathrm{UV}}\gtrsim -22.5$, both semi-analytical models (e.g., Mason et al. 2015a) and cosmological simulation (e.g., Trac et al. 2015; Cowley et al. 2018; Yung et al. 2018).

At the brightest magnitude bin of z ∼ 10, however, the estimated density is ≳1 dex above what theoretical models expect. While our sample size is limited, and follow-up studies are necessary to differentiate those candidates from the contamination, this may highlight factors that were missed in previous models, such as evolution of the shape of LFs. The contribution of AGNs, which possibly boost the observed light of our candidates, would also be worth investigating.

Our magnitude range is limited to the bright end, and thus deriving robust LFs with the present data alone is challenging. Still, it is worth investigating how much our new candidates influence LFs derived in previous studies that focused on faint objects but with limited volume, where bright objects could be easily missed. We attempt here to combine our z ∼ 9 candidates with those in Ishigaki et al. (2018) and investigate the impact on the best-fit parameters.¹⁸

The z ∼ 9 candidates in Ishigaki et al. (2018) are selected by color–color selections from all the Hubble Frontier Fields (∼56 arcmin² in total), similar to our scheme. We take the data points of stepwise LFs in their Figure 4, which already takes into account the effect of lens magnification. While their samples include those with photometric redshifts inconsistent with the selection, we still use their values to reproduce their best-fit parameters of LF. The error, dominated by Poisson error, is recalculated in the same manner as ours.

We first fit only their points with a similar technique, the Markov chain Monte Carlo (MCMC) method, to see whether our fitting method reproduces their Schechter function fit:

$$\begin{eqnarray}\begin{array}{rcl}\phi ({M}_{\mathrm{UV}}) & = & \displaystyle \frac{\mathrm{ln}10}{2.5}{\phi }^{* }\times {10}^{0.4(\alpha +1)({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}\\ & & \times \,\exp [-{10}^{-0.4({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}].\end{array}\end{eqnarray} \tag{ 5 }$$

We use the emcee code (Foreman-Mackey et al. 2013), assuming flat priors in the following ranges: $\alpha \in [-10\,:0]$, ${M}^{* }\in [-30\,:0]$ mag, and $\mathrm{log}{\phi }^{* }\in [-10\,:0]$ Mpc⁻³. The calculation is repeated ${N}_{\mathrm{mc}1}={\rm{10,000}}$ times, and the first half realization is discarded to avoid the bias from initial parameter values. The best-fit values and uncertainties (16/50/84 percentiles) are calculated from the rest of the chain.

When all the parameters are set as free, the fitting results are unconstrained with large uncertainty, as also seen in Ishigaki et al. (2018), because the knee of the LF (${M}_{\mathrm{UV}}\lesssim 21$ mag) is not sufficiently sampled by their data only. We then follow Ishigaki et al. (2018) and fix two of the parameters ($\alpha =-1.96$ and ${M}_{\mathrm{UV}}^{* }=-20.35$), finding consistent values for the others. This ensures that results with our fitting method are comparable.

We then add our data and fit with the Schechter function, but with a few updates in the treatment of uncertainties. We run the MCMC fitting routine introduced above ${N}_{\mathrm{mc}2}=300$ times by fluctuating the data, as described in the following.

First, we fluctuate the observed flux within the random flux error and the systematic uncertainty from the magnification model in each iteration. While the latter uncertainty is relatively small for our candidates, some of those of Ishigaki et al. (2018) have large uncertainties in the magnification (by a factor of ∼3) because of complicated cluster lens modeling.¹⁹

Second, we take account of the redshift uncertainty. In Ishigaki et al. (2018), they did not have a selection criterion with photometric redshifts, while our final sample does include one (Section 2.5). To make the two samples consistent, we fluctuate redshifts within their 1σ uncertainty range in each iteration and select those with $7.7\lt z\lt 9.7$ (95% confidence interval for our z ∼ 9 candidates) when calculating the number density. Absolute magnitudes are also recalculated based on fluctuated redshifts in this step.

Finally, we refine the UV magnitude grid when calculating the number density, because different grids change the stepwise number density and can affect the fitting parameters (Schmidt et al. 2014). We set the magnitude grid as ${M}_{\mathrm{UV}}\,\in [-24+\mathrm{rand}(-0.5,0.5)\,:\,-\,13+\mathrm{rand}(-0.5,\,0.5)]$ with a magnitude bin size ${\rm{\Delta }}{M}_{\mathrm{UV}}=\mathrm{rand}(0.3,1.5)$ in each iteration, where $\mathrm{rand}(a,b)$ is a random float value taken from the range between a and b. The effective volume is interpolated to the refined magnitude grid. Skipping the third process would underestimate the uncertainty down to ∼30%, depending on the size of the magnitude bin.

The result from combining our new observations with the literature data is shown in Figure 6, where we find that all three parameters are constrained within physically meaningful ranges (Table 4). We present the 50/16/84th percentiles taken from the synthesized MCMC chain (${N}_{\mathrm{mc}1}\times {N}_{\mathrm{mc}2}$) as the best-fit parameters. Compared to the fixed values of Ishigaki et al. (2018) ($\alpha =-1.96$ and ${M}_{\mathrm{UV}}^{* }=-20.35$), our value for α is smaller (−2.1), and ${M}_{\mathrm{UV}}^{* }$ is ∼0.6 mag smaller (−21.0), which can be understood from the fact that our data constrain the bright end at a lower number density than previously possibly.

Zoom In Zoom Out Reset image size

Table 4.  Fitting Parameters for Schechter Luminosity Function and Cosmic Star Formation Rate Density at z ∼ 9

$\mathrm{log}{\phi }^{* }$	${M}_{\mathrm{UV}}^{* }$	α	$\mathrm{log}{\rho }_{\mathrm{UV}}$		$\mathrm{log}{\psi }_{* }$
(Mpc−3 mag−1)	(mag)		(erg s−1 Hz−1 Mpc−3)		(${M}_{\odot }$ yr−1 Mpc−3)
			(${M}_{\mathrm{UV}}\lt -17$)	(${M}_{\mathrm{UV}}\lt -15$)	(${M}_{\mathrm{UV}}\lt -17$)	(${M}_{\mathrm{UV}}\lt -15$)
$-{4.24}_{-0.93}^{+0.56}{\ }_{-0.38}^{+0.44}$	$-{21.01}_{-1.35}^{+0.69}{\ }_{-0.50}^{+0.58}$	$-{2.06}_{-0.29}^{+0.31}{\ }_{-0.15}^{+0.27}$	${25.33}_{-0.12}^{+0.01}{\ }_{-0.09}^{+0.09}$	${25.55}_{-0.24}^{+0.10}{\ }_{-0.15}^{+0.15}$	$-{2.61}_{-0.12}^{+0.01}{\ }_{-0.09}^{+0.09}$	$-{2.39}_{-0.24}^{+0.10}{\ }_{-0.15}^{+0.15}$

Note. 50th and 16/84th percentiles are taken from the MCMC realization as the best-fit values and their uncertainties. Associated errors are random photometric errors and systematic errors from the binning size in the UV absolute magnitude. We use ${\kappa }_{\mathrm{UV}}=1.15\times {10}^{-28}\,{M}_{\odot }$ yr⁻¹/(erg s⁻¹ Hz⁻¹) (Madau & Dickinson 2014) to convert from the UV luminosity density (${\rho }_{\mathrm{UV}}$) to the star formation rate density (${\psi }_{* }$).

Download table as:  ASCIITypeset image

With the best-fit Schechter parameters, we estimate the cosmic UV luminosity density,

$$\begin{eqnarray}&&{\rho }_{\mathrm{UV}}={\int }_{-\infty }^{{M}_{\mathrm{lim}.}}{{dM}}_{\mathrm{UV}}L({M}_{\mathrm{UV}})\phi ({M}_{\mathrm{UV}})\end{eqnarray} \tag{ 6 }$$

where $L({M}_{\mathrm{UV}})$ is the UV luminosity at a given UV absolute magnitude. As summarized in Table 4, we find $\mathrm{log}{\rho }_{\mathrm{UV}}\sim 25.3$ (25.6) erg s⁻¹ Hz⁻¹ Mpc⁻³ with ${M}_{\mathrm{lim}}=-17$ (−15), or ${\rm{l}}{\rm{o}}{\rm{g}}{\psi }_{\ast }\sim -2.6$ (−2.4) M_⊙ yr⁻¹ Mpc⁻³ when the conversion in Madau & Dickinson (2014) is applied.

The derived value is consistent with that derived in Ishigaki et al. (2018, $\mathrm{log}{\psi }_{* }\sim -2.7\,{M}_{\odot }$ yr⁻¹ with $\alpha =-1.96$ and ${M}_{\mathrm{UV}}^{* }=-20.35$) or slightly larger than that in Oesch et al. (2013, $\mathrm{log}{\psi }_{* }\sim -2.9$ with $\alpha =-1.73$ and $\mathrm{log}{\phi }^{* }=-2.94$), but in excellent agreement with McLeod et al. (2016) and theoretical models (e.g., Tacchella et al. 2013; Mason et al. 2015b, see also Figure 5 of Oesch et al. 2018). The luminosity density derived here lies between those claimed at lower redshifts (e.g., z ∼ 8 in Bouwens et al. 2015b) and higher redshift (z ∼ 10 in Oesch et al. 2018).

Our updated value of the UV luminosity density does not qualitatively change the conclusion in Ishigaki et al. (2018), where they found a linear relation in z–$\mathrm{log}{\rho }_{\mathrm{UV}}$ up to z ∼ 9 (see also Oesch et al. 2013; Bouwens et al. 2016; McLeod et al. 2016), which may then turn into an accelerated decrease at yet higher redshift (Oesch et al. 2018).

We also fit the number density with a double power law, motivated by recent studies at lower redshift (e.g., Ono et al. 2018):

$$\begin{eqnarray*}\begin{array}{rcl}\phi ({M}_{\mathrm{UV}}) & = & \displaystyle \frac{\mathrm{ln}10}{2.5}{\phi }^{* }\times [{10}^{0.4(\alpha +1)({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}\\ & & +{{10}^{0.4(\beta +1)({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}]}^{-1}.\end{array}\end{eqnarray*}$$

We set the same prior range as above for α, ${M}_{\mathrm{UV}}^{* }$, and ${\phi }^{* }$, and set $\beta \in [-10\,:\alpha ]$. The fit shows unconstrained β and large uncertainties for the other parameters when all the parameters are free. When the characteristic magnitude is fixed (${M}^{* }\,=-20.35$), the other parameters return some constrained values, but still with relatively large uncertainties ($\mathrm{log}{\phi }^{* }=-{3.48}_{-0.11}^{+0.11}$, $\alpha =-{1.92}_{-0.17}^{+0.17}$, and $\beta =-{5.80}_{-1.89}^{+2.00}$). Thus, we conclude that we are still unable to rule out this functional form, and that additional large-area surveys would be highly beneficial.