# Bright $z\sim 9$ Galaxies in Parallel: The Bright End of the Rest-UV
Luminosity Function from HST Parallel Programs

Micaela B. Bagley Department of Astronomy, The University of Texas at Austin,
Austin, TX, USA  Steven L. Finkelstein Department of Astronomy, The University
of Texas at Austin, Austin, TX, USA  Sofía Rojas-Ruiz Fellow of the
International Max Planck Research School for Astronomy and Cosmic Physics at
the University of Heidelberg (IMPRS–HD)  Max-Planck-Institut für Astronomie,
Königstuhl 17, D-69117, Heidelberg, Germany  James Diekmann  Department of
Astronomy, The University of Texas at Austin, Austin, TX, USA  Keely D.
Finkelstein Department of Astronomy, The University of Texas at Austin,
Austin, TX, USA  Mimi Song University of Massachusetts, Amherst, MA, USA
Casey Papovich George P. and Cynthia Woods Mitchell Institute for Fundamental
Physics and Astronomy, Department of Physics and Astronomy, Texas A&M
University, College Station, TX, USA  Rachel S. Somerville  Center for
Computational Astrophysics, Flatiron Institute, NY, USA  Ivano Baronchelli
INAF- Istituto di Radioastronomia - Italian ALMA Regional Centre, via Gobetti
101, 40129 Bologna, Italy  Y.Sophia Dai Chinese Academy of Sciences South
America Center for Astronomy (CASSACA), National Astronomical Observatories of
China (NAOC), 20A Datun Road, Beijing, 100012, China

###### Abstract

The abundance of bright galaxies at $z>8$ can provide key constraints on
models of galaxy formation and evolution, as the predicted abundance varies
greatly when different physical prescriptions for gas cooling and star
formation are implemented. We present the results of a search for bright
$z\sim 9-10$ galaxies selected from pure-parallel Hubble Space Telescope
imaging programs. We include 132 fields observed as part of the Brightest of
Reionizing Galaxies survey, the Hubble Infrared Pure Parallel Imaging
Extragalactic Survey, and the WFC3 Infrared Spectroscopic Parallel survey.
These observations cover a total of 620 arcmin2, about 70% of which is also
covered with Spitzer Space Telescope infrared imaging. We identify thirteen
candidate galaxies in the range $8.3<z<11$ with $24.5<m_{H}<26.5$
($-22.9<M_{UV}<-21.2$). This sample capitalizes on the uncorrelated nature of
pure parallel observations to overcome cosmic variance and leverages a full
multi-wavelength selection process to minimize contamination without
sacrificing completeness. We perform detailed completeness and contamination
analyses, and present measurements of the bright end of the UV luminosity
function using a pseudo-binning technique. We find a number density consistent
with results from Finkelstein et al. (2022) and other searches in HST parallel
fields. These bright candidates likely reside in overdensities, potentially
representing some of the earliest sites of cosmic reionization. These new
candidates are excellent targets for follow-up with JWST, and four of them
will be observed with the NIRSpec prism in Cycle 1.

††facilities: HST (WFC3), Spitzer (IRAC), MAST, IRSA††software: Astropy
(Astropy Collaboration et al., 2013), Photutils (v0.6; Bradley et al., 2019),
SciPy (Jones et al., 2001–), Scikit-learn (Pedregosa et al., 2011), Source
Extractor (Bertin & Arnouts, 1996), Galfit (Peng et al., 2010), EAZY (Brammer
et al., 2008)

##  1 Introduction

The study of galaxies in the very early universe, and particularly of UV-
bright galaxies, provides key input to models of galaxy formation. The number
density of UV-bright galaxies is affected by factors including feedback
processes (e.g., Somerville et al., 2008; Bower et al., 2012), dust
attenuation (e.g., Finkelstein et al., 2012a; Vogelsberger et al., 2020), the
build-up of dark matter halos, and star formation efficiency. Constraining the
abundance of UV-luminous galaxies thus directly constrains the fundamental
physics of star formation in the early Universe. By probing ever-larger
volumes, recent ground and space-based surveys are providing a progressively
more complete analysis of the abundance of bright galaxies. These studies are
showing that the characteristic luminosity ($L^{*}$) of the rest-frame UV
luminosity function only shallowly evolves (if at all) in the range $z=4-8$
(Bowler et al., 2014, 2015, 2020; Bouwens et al., 2015; Finkelstein et al.,
2015b). As the UV luminosity function probes recent star formation, these
results imply that star-formation rates evolve more slowly than the halo mass
function, which would indicate that galaxies are efficient at converting gas
into stars even out to these redshifts (e.g., Behroozi & Silk, 2015;
Finkelstein et al., 2015c; Yung et al., 2019a, b, though see Stefanon et al.
2021). If this observational trend continues to higher redshifts, it will
provide a significant challenge to current galaxy formation models.

Considerable effort has been devoted to modeling galaxy evolution at early
cosmic times ($z\gtrsim 6$), via methods such as hydrodynamical simulations
(e.g., Gnedin, 2016; Wilkins et al., 2017), (semi-)analytic models (e.g.,
Somerville et al., 2015; Yung et al., 2019a, b), and empirical models (e.g.,
Mason et al., 2015a; Behroozi et al., 2020). Yet predictions for the density
of galaxies at $z>9$ can differ dramatically. As the predicted number of UV-
bright galaxies is sensitive to the assumed relationship between molecular-gas
density and the star-formation-rate density in cosmological simulations,
robust observations in this epoch will provide powerful constraints on these
physical processes. Analyses at $z>8$ continue to push the limits of what is
possible with the current data. The majority of candidates at these redshifts
(e.g., Coe et al., 2013; Oesch et al., 2013, 2014; Bouwens et al., 2015, 2016,
2019; Finkelstein, 2016) were selected from a small set of well-studied fields
such as the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey
(CANDELS; Grogin et al., 2011; Koekemoer et al., 2011), the Hubble Ultra Deep
Field (HUDF; Beckwith et al., 2006), and the Hubble Frontier Fields (Lotz et
al., 2017). Results measured from a few independent pointings are susceptible
to cosmic variance, which can dominate the uncertainties at the bright end of
the UV luminosity function. Thus, the density of UV-bright galaxies is
presently ill-constrained at $z>8$, despite its importance for understanding
high-redshift star formation and galaxy evolution.

Pure parallel programs – where random field pointings result in completely
independent, uncorrelated observations – provide a key opportunity to address
these outstanding questions. Often shallower than targeted surveys, pure
parallel programs are well-suited to identifying bright high-redshift
galaxies: a population of rare sources that may not be fully sampled by
surveys covering the same area in contiguous fields. Parallel observing can
therefore provide better constraints on the bright end of the UV luminosity
function and is complementary to deeper surveys that probe the fainter
population.

Pure parallel programs such as the Brightest of Reionizing Galaxies (BoRG, PI:
Trenti; General Observer [GO] programs 11700, 12572) and the Hubble Infrared
Pure Parallel Imaging Extragalactic Survey (HIPPIES, PI: Yan; GO 11702, 12286)
were indeed designed for just this purpose. Using WFC3 imaging in at least
four filters to probe the Lyman break at $z\sim 8$, these surveys drastically
increased the known sample of bright galaxies and provided crucial constraints
on the luminosity function (Trenti et al., 2011, 2012; Bradley et al., 2012;
Schmidt et al., 2014a). The independent pointings significantly reduce the
effect of cosmic variance from $\gtrsim$20% to $<$1% (Trenti & Stiavelli,
2008; Bradley et al., 2012), providing comparable constraining power as a
similar-depth survey with twice the area but in a single pointing (Bradley et
al., 2012). Studies using the BoRG[z910] fields, which included
$JH_{140}$-band imaging, find a relatively high number density of bright
$z\sim 9-10$ galaxies (Calvi et al., 2016; Morishita et al., 2018; Rojas-Ruiz
et al., 2020), potentially indicating more efficient star formation at these
early times (similar results have been seen in the CANDELS fields; Finkelstein
et al., 2022). Yet significant scatter also remains, demonstrating the
systematic uncertainty resulting from different sample selection techniques
(i.e., strict color cuts versus photometric redshift selection). While JWST
will quickly shed light on this tension, it is not optimized to conduct wide-
field surveys for new bright, high-redshift galaxies. The power of JWST will
lie in targeted follow-up for redshift confirmation and source
characterization.

To that end, we have performed a new independent selection for $z=9-10$
galaxies in a subset of available HST pure parallel fields, improving on
previous efforts by including all available HST and Spitzer/IRAC photometry in
our photometric-redshift selection. We also included in our search 45 fields
($\sim$160 arcmin2) from the WFC3 Infrared Spectroscopic Parallel survey
(WISP, PI: Malkan; GO 12283, 12902, 13352, 13517, 14178), the first time these
fields have been searched for $z\sim 9$ candidates. While many have performed
searches for $z>9$ galaxies in the BoRG/HIPPIES fields (Bernard et al., 2016;
Calvi et al., 2016; Morishita et al., 2018; Bridge et al., 2019; Rojas-Ruiz et
al., 2020; Morishita et al., 2020; Morishita, 2021; Roberts-Borsani et al.,
2021), the results at the bright end of the UV luminosity function remain
uncertain due in part to varied selection techniques probing close to the
image detection limits. In fact, these teams find different samples of galaxy
candidates _in the same fields_ , highlighting how sensitive HST results at
$z\gtrsim 9$ are to photometric measurements, noise characterizations, and
sample selections. In our search, we employ a full multi-wavelength selection
(independent of a simple $J_{125}$$-$$H_{160}$ color cut) including
Spitzer/IRAC photometry where available, with detailed noise calculations and
source vetting to remove contaminants – resulting in a larger, more complete
sample. This sample includes some of the brightest candidates at these
redshifts, potentially representing some of the most massive galaxies to form
$\lesssim$500 Myr after the Big Bang.

This paper is organized as follows. In Section 2, we present the pure parallel
imaging datasets and describe our photometric measurements. We describe our
sample selection criteria in Section 3, including a detailed discussion of how
we vetted each candidate, performing a thorough check for cases of persistence
(Section 3.3) and a careful visual inspection (Section 3.4). We incorporate
available Spitzer/IRAC imaging of each candidate in Section 3.5, and check for
stellar contamination in Section 3.7. In Section 4, we present our 13 high-
redshift candidate galaxies and compare them to previous studies in the same
HST fields. We address the possibility of low-redshift contamination in
Section 5, and use our sample to infer the density of galaxies at the bright
end of the UV luminosity function in Section 6. Finally, we discuss
implications for our sample of bright, high-redshift galaxies in Section 7 and
summarize in Section 8. Throughout this paper we assume a $\Lambda$CDM
cosmology with $\Omega_{M}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km s-1
Mpc-1 and express all magnitudes in the AB system (Oke & Gunn, 1983) unless
otherwise noted.

##  2 Observations

We include observations from multiple programs, all obtained as part of
parallel imaging programs with WFC3111www.stsci.edu/hst/wfc3 (Kimble et al.,
2008) while another HST instrument was in use. Typically these parallel
observations are taken while either the Cosmic Origins Spectrograph (COS;
Froning & Green, 2009) or the Space Telescope Imager and Spectrograph (STIS;
Kimble et al., 1998) are engaged in long integrations of a primary target.
Without control of the telescope pointing, pure parallel programs
(unaffiliated with the primary programs) can accrue observations of
independent and uncorrelated fields during these long primary integrations. In
this paper, we combine data from a subset of the Brightest of Reionizing
Galaxies survey (BoRG; Trenti et al., 2011), the Hubble Infrared Pure Parallel
Imaging Extragalactic Survey (HIPPIES; Yan et al., 2011), the WFC3 Infrared
Spectroscopic Parallel survey (WISP; Atek et al., 2010), and additional
coordinated parallels from the COS Guaranteed Time Observer (GTO) program.

The programs include imaging with both the UVIS and IR cameras on WFC3, as
well as some optical imaging with the Advanced Camera for Surveys
(ACS222www.stsci.edu/hst/acs; Clampin et al., 2000) for HIPPIES observations
from Cycle 18. All fields are observed with the F160W filter, which we use for
detection, and three additional filters covering wavelengths from
$\sim$$0.6-1.4\micron$. True $z\gtrsim 9$ galaxies will be detected in F160W
and at most one additional filter (either F125W or F110W depending on the
survey in question), and the optical imaging is crucial in detecting the
Lyman-$\alpha$ (Ly$\alpha$) break at 1216Å. By $z$$\sim$10 F160W is the only
filter available in these programs that is redward of the Lyman break. The
WFC3/IR camera has a field of view ($123\arcsec\times 134\arcsec$) that is
smaller than both of the optical cameras (WFC3/UVIS: $162\arcsec\times
162\arcsec$, and ACS/WFC: $202\arcsec\times 202\arcsec$) and therefore
determines the total area covered in each field. The native pixel scales of
the cameras are $0\farcs 13$/pixel (IR), $0\farcs 04$/pixel (UVIS), and
$0\farcs 05$/pixel (WFC). The filter coverage and depth achieved in each field
depends on the survey strategies of each program and the length of each
parallel opportunity. The fields we assemble here therefore have inconsistent
depths and filter coverages, which we handle with careful field-specific
simulations to assess the completeness and effective volume probed in each
field.

Table 1: Dataset Summary  Dataset | $N_{\mathrm{fields}}$ | $A_{\mathrm{total}}$ | Pixel scale
---|---|---|---
|  | (arcmin2) | (arcsec/pixel)
BoRG[$z$8] | 64 | 308.8 | 0.08
HIPPIES (GO 12286) | 23 | 108.3 | 0.10
WISP | 45 | 203.9 | 0.13
Figure 1: The positions of all 132 fields we consider in our analysis. Purple
circles indicate BoRG[$z$8] fields, which include 34 fields from the BoRG
survey (PI: Trenti), 8 HIPPIES fields (PI: Yan,GO 11702), and 22 parallel
fields observed as part of the COS GTO program (PI: Green). The 23 HIPPIES
fields from GO 12286 are shown as orange squares, and the 45 WISP fields are
shown as blue triangles. (Symbol sizes are not to scale.) The independent,
uncorrelated nature of the fields – ideal for studies of rare, bright galaxies
– is evident from their locations shown here.

In total, we consider 132 WFC3 parallel fields covering $\sim$620 arcmin2: 64
BoRG[$z$8] fields (comprised of 34 fields from the BoRG survey, 8 HIPPIES
fields from GO program 11702, and 22 parallel pointings from the COS GTO
program), 23 additional HIPPIES fields (GO 12286), and 45 fields from the WISP
survey. The locations of all 132 fields are shown in Figure 1, where the
independent nature of the pointings is clear. We provide a summary of these
datasets in Table 1 and field-specific information about all 132 fields,
including filter coverage and imaging depths, in Table 2. We measure the area
covered in each pointing by summing the pixels in the F160W weight maps with
values greater than 100 counts/sec ($=\mathrm{RMS}<0.1$), as this allows us to
account for the area lost to cosmic rays, detector artifacts, or other bad
pixels. The imaging depths presented in Table 2 are the $5\sigma$ limiting
magnitudes as measured in circular apertures with radius $r=0\farcs 2$ via the
following method. We place $\sim$5000$-$10000 uncorrelated apertures randomly
across the background of each image, avoiding source flux and bad pixels. The
limiting magnitude quoted in each filter is found by measuring the standard
deviation of the distribution of fluxes from all apertures, and thus
represents an average depth across the full image. We note that we have not
performed an aperture correction to account for the fractional flux enclosed
by the $r=0\farcs 2$ apertures. The measured depth can vary across a field by
up to $\pm 0.3$ magnitudes from the median, and so we rely on the noise
measured locally around each source (see Section 2.4.1) to determine detection
significances that we use when selecting high-redshift candidate galaxies (see
Section 3.2).

We briefly explain the characteristics and data reduction of each dataset in
the following sections. We hereafter refer to the HST filters using a letter
that specifies the bandpass followed by the three-number identification: i.e.,
$I_{600}$ for the long pass filter F600LP, $H_{160}$ for F160W.

Table 2: Parallel Field Coverage and Depths

Field | RA | Dec |  $V_{606}$aaThe 5$\sigma$ limiting magnitudes are measured for each filter in circular apertures with $r=0\farcs 2$. | $I_{600}$ | $I_{814}$ | $Y_{098}$ | $Y_{105}$ | $J_{110}$ | $J_{125}$ | $H_{160}$ | PI | Program | AreabbThe area in arcmin2 is measured as $H_{160}$ pixels with weights of $>100$ (RMS $<0.1$).
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) | (mag) |  |  | (arcmin2)
Par0110$-$0224accFor a few of the BoRG[$z$8] and HIPPIES fields, two overlapping HST pointings have been combined to form a single field, resulting in varied filter coverage and depths. In these cases, we have separated the field into two components, a and b, and consider them separately in our analysis. See the text for details.  | 01:10:04.4 |  $-$02:25:04.6 | 27.54 | 27.09 |  $\cdots$ | 27.99 |  $\cdots$ |  $\cdots$ | 27.74 | 27.42 | Trenti | 11700 | 6.45
Par0110$-$0224b | 01:10:04.4 |  $-$02:25:04.6 | 27.51 |  $\cdots$ |  $\cdots$ | 27.60 |  $\cdots$ |  $\cdots$ | 27.55 | 27.18 | Trenti | 11700 | 8.22
Par0214+1255 | 02:13:37.8 |  $+$12:54:51.6 |  $\cdots$ | 26.09 |  $\cdots$ | 26.70 |  $\cdots$ |  $\cdots$ | 26.64 | 26.49 | Yan | 11702 | 4.62
Par0228$-$4102 | 02:27:56.0 |  $-$41:01:34.8 |  $\cdots$ | 27.16 |  $\cdots$ | 27.52 |  $\cdots$ |  $\cdots$ | 27.39 | 27.16 | Green | 11541 | 4.64
Par0240$-$1857 | 02:40:26.9 |  $-$18:57:19.6 |  $\cdots$ | 27.05 |  $\cdots$ | 27.51 |  $\cdots$ |  $\cdots$ | 27.38 | 27.13 | Green | 11541 | 4.63
Par0005+1607 | 00:05:39.5 |  $+$16:07:05.3 |  $\cdots$ | 27.06 |  $\cdots$ |  $\cdots$ | 27.17 |  $\cdots$ | 27.47 | 27.05 | Yan | 12286 | 4.62
Par0241+0715 | 02:41:48.3 |  $+$07:15:35.5 |  $\cdots$ | 26.76 |  $\cdots$ |  $\cdots$ | 27.16 |  $\cdots$ | 27.21 | 26.88 | Yan | 12286 | 4.62
Par0245+1051 | 02:45:25.7 |  $+$10:51:57.9 |  $\cdots$ | 26.80 |  $\cdots$ |  $\cdots$ | 27.34 |  $\cdots$ | 27.29 | 26.91 | Yan | 12286 | 4.62
Par0259+0032 | 02:59:33.4 |  $+$00:32:13.4 |  $\cdots$ | 26.71 |  $\cdots$ |  $\cdots$ | 27.08 |  $\cdots$ | 27.11 | 26.50 | Yan | 12286 | 4.60
Par0713+7405 | 07:13:22.2 |  $+$74:06:03.3 |  $\cdots$ | 27.06 |  $\cdots$ |  $\cdots$ | 27.03 |  $\cdots$ | 27.33 | 27.18 | Yan | 12286 | 2.89
Par96 | 02:09:24.3 |  $-$04:43:50.2 | 27.48 |  $\cdots$ | 26.88 |  $\cdots$ |  $\cdots$ | 27.94 |  $\cdots$ | 26.97 | Malkan | 12283 | 4.59
Par104 | 10:05:25.1 |  $+$01:29:43.5 | 27.02 |  $\cdots$ | 26.44 |  $\cdots$ |  $\cdots$ | 27.78 |  $\cdots$ | 26.73 | Malkan | 12283 | 4.58
Par257 | 02:46:37.0 |  $-$30:32:05.5 | 26.90 |  $\cdots$ | 26.47 |  $\cdots$ |  $\cdots$ | 27.51 |  $\cdots$ | 26.42 | Malkan | 12902 | 4.53
Par296 | 10:27:11.4 |  $+$18:35:31.1 | 27.04 |  $\cdots$ | 26.63 |  $\cdots$ |  $\cdots$ | 27.19 |  $\cdots$ | 26.40 | Malkan | 12902 | 4.54
Par298 | 09:21:30.5 |  $+$45:07:9.5 | 26.97 |  $\cdots$ | 26.44 |  $\cdots$ |  $\cdots$ | 26.99 |  $\cdots$ | 26.24 | Malkan | 12902 | 4.26

Note. — Five representative fields from each dataset are presented here. This
table is available in its entirety in machine-readable format in the
electronic version of the paper.

Figure 2: The filter coverage of each dataset, showing that of BoRG[$z$8] (top
panel), HIPPIES GO 12286 (middle), and WISP (bottom). The $V_{606}$ and
$I_{600}$ profiles in the top row have dashed outlines to indicate that the
BoRG[$z$8] fields are observed with one or the other of these filters. In
place of $J_{125}$, the WISP survey observed with the broader $J_{110}$,
making WISP fields less sensitive to galaxies at $z\sim 9$ than at $z\sim 10$.
For reference, we show a model galaxy spectrum at $z=9$ and $z=10$ with a weak
Ly$\alpha$ emission line. At $z=9$, the Lyman break is near the center of the
$J_{125}$, while by $z>10$ it has redshifted almost completely out of the
filter, resulting in single-band $H_{160}$ detections.

###  2.1 BoRG[$z$8]

We include in our analysis data from the first set of fields observed as part
of the BoRG survey (PI: Trenti; Trenti et al., 2011). These fields were
observed in HST Cycles 17 and 19 (GO 11700, 12572) with the $Y_{098}$,
$J_{125}$, and $H_{160}$ IR filters and the $V_{606}$ UVIS filter (see top row
of Figure 2). The BoRG data releases also include data from programs with
similar observing strategies and filter coverages – HIPPIES (PI: Yan, GO
11702; Yan et al., 2011) and coordinated parallels from the COS GTO team (PI:
Green, GO 11528, 11530, 11534, 11541, 12024, 12025) – though the HIPPIES
fields as well as some of the COS parallels were observed with the $I_{600}$
filter instead of $V_{606}$. Collectively, these three programs cover 71
fields, with data and results described in detail in Trenti et al. (2011),
Bradley et al. (2012), and Schmidt et al. (2014a). These surveys were
optimized to detect $z$$\sim$8 galaxies as $Y_{098}$-dropouts with detections
in both $J_{125}$ and $H_{160}$ and are therefore referred to collectively as
BoRG[$z$8]. This designation distinguishes the first set of BoRG observations
with those from later HST cycles that included the $JH_{140}$ filter. Imaging
with $JH_{140}$ allows for the detection of $z\gtrsim 9$ galaxies as
$J_{125}$-dropouts with two-band detections redward of the break. These
BoRG[$z$910] fields have been used to identify $z\gtrsim 9$ candidates in
multiple studies (Calvi et al., 2016; Morishita et al., 2018, 2020; Morishita,
2021; Bridge et al., 2019; Rojas-Ruiz et al., 2020; Roberts-Borsani et al.,
2021). In this paper we focus on the BoRG[$z$8] fields, and perform an
analysis similar to that presented in Rojas-Ruiz et al. (2020).

We use the reduced images and weight maps available from the BoRG data
Delivery 3 from MAST333https://archive.stsci.edu/prepds/borg. These images are
all pixel-aligned and on a scale of 0$\farcs$08/pixel. We convert the weight
maps (actually inverse variance maps output by the WFC3 pipeline, CALWFC3) to
RMS maps for use with source detection via
$\mathrm{RMS}=1/\sqrt{\mathrm{WHT}}$. Pixels with zero weight are assigned a
value of 104 in the RMS maps.

After visually inspecting all 71 BoRG[$z$8] fields, we chose to remove seven
fields from consideration in our analysis. Three of the fields – Par0540-6409,
Par0553-6405, and Par1815-3244 – appear to cover star clusters, resulting in a
very high density of point sources across the full image. Two others –
Par1014-0423 and Par1230+0750 – are heavily contaminated by strong detector
persistence, likely resulting from observations of bright star clusters
immediately preceeding the BoRG[$z$8] exposures. We remove the remaining two –
Par0751+2917 and Par1209+4543 – because they were also observed in Cycle 22
with additional filters (F350LP and $JH_{140}$) and so are part of the
BoRG[$z$910] survey. One field, Par0110-0224, was covered by both the BoRG and
HIPPIES GO 11702 programs and is therefore partially covered by five filters
(including both $V_{606}$ and $I_{600}$). We separated this field into two
parts – one with and one without $I_{600}$ – which are labeled $a$ and $b$,
respectively, in Table 2. We count this as a single field in the total (of the
132), but treat them as separate fields in all other aspects of our analysis.
Specifically, we perform separate completeness simulations to account for the
varying filter coverage and depth across the field. Finally, a few fields are
partially affected by significant detector issues or other contamination,
including scattered light from Earth’s limb, highly variable background across
part of the images, and extremely strong Dragon’s Breath (Fowler et al., 2017)
in the UVIS filters. In these cases, we mask out the affected regions in all
filters and weight maps, therefore accounting for the corresponding loss in
total survey area. The field areas listed in Table 2 reflect the unaffected,
unmasked portions of these fields. In total, we include 64 BoRG[$z$8] fields
in our analysis, covering $\sim$308 arcmin2.

###  2.2 HIPPIES

The HIPPIES observations from Cycle 17 were reduced and released as part of
BoRG[$z$8] (see previous section), due to the similar filter sets and
observing strategies of the two programs. However, HIPPIES observations from
Cycle 18 (PI: Yan, GO 12286) were not included, and so we count them here as a
separate dataset. This program used the near-IR filter $Y_{105}$ rather than
the medium band $Y_{098}$, and also used $I_{600}$ as the bluest filter.

We downloaded the FLTs and supporting files for this dataset from MAST and
created reduced mosaics following the same procedure outlined in Rojas-Ruiz et
al. (2020), which we briefly summarize here. The final, drizzled images are
produced using a pipeline custom built for use with HIPPIES data, with the
modifications necessary to account for the added challenges of undithered,
pure parallel observations. The pipeline uses an MCMC sampler to align
individual images based on detected source positions. It then creates mosaic
images with a pixel scale of 0$\farcs$1/pixel using the MultiDrizzle software
package (Koekemoer et al., 2003). The output weight maps are then scaled by
the average amplitude of the correlation measured in blank regions of the
science images in order to produce RMS maps.

We rejected a handful of pointings that are severely impacted by features such
as scattered Earth light, and we masked portions of two others. We also
separated field Par2346-0021 into two parts due to inconsistent filter
coverage and imaging depth across the field, similar to that described for
Par0110-0224 for the BoRG[$z$8] dataset. We include 23 fields ($\sim$108
arcmin2) from the HIPPIES GO program 12286.

###  2.3 WISP

Similar to the programs described above, WISP (PI Malkan; Atek et al., 2010)
obtained imaging of uncorrelated fields through pure parallel observations.
However, WISP also acquired slitless spectroscopy using WFC3’s two near-
infrared grisms: G102 and G141. Primarily a spectroscopic program, WISP
devoted the majority of available integration time to grism observations and
obtained imaging mainly to aid in extracting and calibrating spectra from the
slitless grism images. Yet for long parallel opportunities consisting of
$\geq$5 orbits, WISP pointings were observed with four imaging filters along
with both grisms. We include 45 of these “deep” WISP pointings ($\sim$204
arcmin2) in our analysis, representing the first search for $z>8$ galaxies in
WISP observations. The WISP pointings are observed with two IR filters
($J_{110}$ and $H_{160}$) and the UVIS filter $I_{814}$. Most of the fields
are also observed with the UVIS filter $V_{606}$.

The main difference between the WISP observations we include here and those of
the other datasets is the use of the $J_{110}$ filter rather than $J_{125}$.
The considerable width of the $J_{110}$ filter (FWHM$\sim$5000Å, compared to
$\sim$3000Å for $J_{125}$) makes it challenging to identify sources at $z\sim
8-9$ via the dropout technique. Galaxies at redshifts of $z\sim 7-10$ will all
be detected in $J_{110}$, and the resulting photometric redshift probability
distributions are too broad to reliably distinguish galaxies at the high-
redshift end of this interval from those at the lower end. We account for this
selection effect with our field-specific completeness simulations.
Additionally, all 45 WISP fields were observed at 3.6µm with Spitzer/IRAC,
which helps constrain the photometric redshifts (see Section 3.5). The
detection power for WISP fields returns at $z\sim 10$, when sources drop out
of the $J_{110}$ filter and become single-band $H_{160}$ detections.

The WISP data reduction pipeline is presented by Atek et al. (2010) and
Battisti et al. (in preparation), and uses the WFC3 pipeline CALWF3 with
modifications required to address the specific challenges of parallel
observing. We use the reduced data available through the MAST
archive444archive.stsci.edu/prepds/wisp/ that were produced with version 6.2
of the WISP pipeline. Pixel-aligned images in all filters are available on
$0\farcs 04$ and $0\farcs 13$ pixel scales, and we adopt the $0\farcs 13$
scale for this analysis because it better samples the point spread function –
with the undithered parallel data there are not enough exposures to cover the
IR pixels at the finer resolution.

###  2.4 Source Extraction and Photometry

We detect sources and perform photometry in all three datasets using Source
Extractor (v2.5; Bertin & Arnouts, 1996) with $H_{160}$ as the detection
filter. Our source detection parameters are the result of optimization using
5000 simulated sources inserted randomly into the images from 21 fields. For
these tests, we chose fields that cover the range of depths, filter coverages,
and pixel scales present across all three datasets. For each field, we created
10 realizations of the $H_{160}$ images, adding 200 synthetic sources to each
realization. Half of the synthetic sources were Sérsic profiles with half-
light radii ($r_{1/2}$) and Sérsic indices ($n$) pulled randomly from uniform
distributions in the ranges $0\farcs 05<r_{1/2}<0\farcs 15$ and $0.5<n<4.5$,
respectively. The other half were modeled as two-dimensional Gaussian profiles
with $\sigma_{x}=\sigma_{y}$ pulled uniformly from the range $0\farcs
05<\sigma<0\farcs 1$. As we were looking to optimize detection parameters, the
sources were placed randomly in clean areas of the images, i.e., avoiding bad
pixels and source flux.

We then ran Source Extractor multiple times on each realization, covering a
grid of source detection parameter values. Specifically, we varied the
threshold for detection over the range $0.3\sigma<$ DETECT_THRESH $<3\sigma$
in steps of $0.1\sigma$ and the minimum number of pixels above the threshold
required for detection over the range $3<$ DETECT_MINAREA $<10$ pixels. We
also inverted each image realization and ran the same grid of parameter values
to quantify the number of negative pixel groupings that would be detected as
real sources. In this way, we characterized the contamination rates of
spurious sources that would be detected as contiguous clusters of noise for
each set of detection parameters. We explored the fraction of recovered
simulated sources as well as the fraction of spurious detections as a function
of $H_{160}$ magnitude, choosing the parameter values that maximize the
recovered fraction ($\geq$80%) while minimizing the spurious fraction ($<$5%)
down to $m_{H160}=27$: DETECT_THRESH = 1.1$\sigma$ and DETECT_MINAREA = 4
pixels. Finally, we found that the same set of detection parameters would
optimize these fractions for all three datasets for $m_{H160}<27$. Any
differences in the recovery and contamination rates for the datasets occur at
fainter magnitudes than we will consider in this paper. Therefore, although
the images are on different pixel scales, we use one consistent set of Source
Extractor detection parameters for all fields.

For each field, we create catalogs using Source Extractor in dual image mode
with $H_{160}$ as the detection filter and the RMS maps for the weights used
in source detection and analysis. Photometry is performed in both circular and
elliptical apertures, with each measurement serving a different purpose in our
sample selection (see Section 3.2). We use small circular apertures with
radius $r=0\farcs 2$ to measure the detection significance of sources in all
filters, because they capture the light in the core of sources while
minimizing the contribution of off-source noise. We use small elliptical Kron
(1980)-like apertures to measure source colors, setting the Source Extractor
parameter PHOT_AUTOPARAMS to (1.2,1.7), where the first number is the scaling
parameter for the ellipse and the second is the minimum radius below which
elliptical apertures are converted to circular. As shown by Finkelstein et al.
(2010), these choices are motivated by the need to maximize both the recovered
flux and the signal-to-noise (S/N) for faint, compact galaxies. Additionally,
the smaller aperture scaling parameter (here 1.2 compared to the Source
Extractor default 2.5) reduces the number of small sources with elliptical
apertures that are artifically stretched to larger sizes by bright neighbors
in crowded regions. In order to similarly determine the optimal parameter
choices for these parallel imaging datasets, we tested the recovered S/N of
the simulated sources inserted into the 21 fields described above over a grid
of PHOT_AUTOPARAMS. The recovered S/N of sources in the $H_{160}$ images peaks
for parameter values in the range ($1-1.5$,$1.3-2$) with no clear winner, and
so we adopt (1.2,1.7), following Bouwens et al. (2010, 2021) and Finkelstein
et al. (2010, 2015b, 2022). We derive aperture corrections for these small
elliptical apertures in Section 2.4.2.

####  2.4.1 Photometric Errors

As the noise in the pixels of drizzled images is partially correlated, the
total noise in an aperture is larger than the sum in quadrature of the noise
in individual (uncorrelated) pixels (Casertano et al., 2000). The flux
uncertainties calculated by Source Extractor are therefore underestimated, and
the S/N measurements based on these flux uncertainties will be artifically
high. We therefore perform our own calculation of the noise in the images
following the process described by Papovich et al. (2016, see also ), which we
briefly summarize here.

Figure 3: An example of the noise calculated for BoRG[$z$8] field Par0846+7654
as a function of aperture size. The noise, $\sigma_{N}$, is calculated as the
normalized median absolute deviation of background flux measured in
uncorrelated circular apertures placed on blank regions of the field. We show
$\sigma_{N}$ for the filters in this field as blue triangles ($I_{600}$),
green diamonds ($Y_{098}$), orange squares ($J_{125}$), and red circles
($H_{160}$). The functional fits to these measurements (equation 1) are
plotted as curves of the corresponding colors. The black dashed and dot-dashed
lines indicate the cases of completely correlated pixels ($\beta=1$) and
completely uncorrelated pixels ($\beta=0.5$), respectively. The noise in these
images is partially correlated in all filters. Neglecting to correct for this
correlation would result in underestimated flux uncertainties and
overestimated S/N’s in our catalog.

In an image with completely uncorrelated pixels, the noise $\sigma_{N}$
measured in an aperture with $N$ pixels is expected to scale as
$\sigma_{N}=\sigma_{b}\times\sqrt{N}$, where $\sigma_{b}$ is the standard
deviation of the pixels containing only the sky background. At the other
extreme, the noise will scale as $\sigma_{N}=\sigma_{b}\times N$ for
completely correlated pixels (Quadri et al., 2007). We thus expect the noise
to scale as $\sigma_{N}\propto N^{\beta}$, with $0.5<\beta<1$ to reflect the
two limiting cases. For our analysis of the noise as a function of aperture
size, we place $\sim 5000$ non-overlapping circular apertures randomly across
each image, using the segmentation and RMS maps to avoid both detected source
flux and bad pixels. The apertures have radii ranging from 0.5 pixels to 12
pixels. We use Photutils (v0.6; Bradley et al., 2019) to measure the flux in
these apertures and estimate $\sigma_{N}$ via the normalized median absolute
deviation (Beers et al., 1990) as a robust estimator of the standard deviation
in the distribution of aperture fluxes. We then fit the following
parameterized function to the relation between $\sigma_{N}$ and $N$:

$\sigma_{N}=\sigma_{b}(\alpha N^{\beta}),$ (1)

where $\alpha$ and $\beta$ are free parameters. For $\sigma_{b}$, we again use
the normalized median absolute deviation to estimate the pixel-to-pixel
variation in background pixels in the image. We show the measured $\sigma_{N}$
and fits for an example field in Figure 3. For each source, filter, and
aperture type, we use equation 1 to calculate a flux error determined by the
number of pixels in the apertures. Finally, we scale the fit by the value of
the RMS map at the position of the source divided by the median of the RMS
map. This scaling maintains any structure present in the RMS map that is
indicative of location-dependent quality differences. We include Poisson
photon errors in this scaling parameter, though we note that including this
shot noise increases the flux uncertainties by only $\sim 3-6$%.

####  2.4.2 Photometric Corrections

We apply three correction factors to the catalog in order to account for (1)
the missing wings of the PSF in the small Kron apertures (i.e., aperture
corrections), (2) any unrecovered source flux due to the source extraction
parameters, and (3) Galactic extinction. These corrections are applied in all
filters to both the measured fluxes and flux uncertainties, thus preserving
the optimized S/N in the smaller apertures while correcting for any missing or
extinguished flux. First, we derive aperture corrections in $H_{160}$ as the
ratio of the fluxes calculated in the small Kron apertures with those using
larger Kron apertures corresponding to PHOT_APERTURES = (2.5,3.5). The median
flux ratio in the full catalog is $\sim$0.69 (i.e., source fluxes are
$\sim$31% higher when measured in the larger apertures), though the aperture
corrections are calculated and applied on a source-by-source basis. Next, we
use the simulated sources from our completeness simulations (see Section 6.2)
to quantify any unrecovered source flux. We find that the median recovered
$H_{160}$ brightness of sources detected with S/N$\geq$5 in all datasets is
$\sim$0.14 magnitudes fainter than the input brightness after application of
the aperture correction. We apply this correction factor to the total
(aperture-corrected) fluxes and uncertainties, obtaining a final estimate of
the total source flux. Lastly, we correct for Galactic extinction using the
color excess $E(B-V)$ for each parallel field taken from the IRSA Galactic
Dust Reddening and Extinction
calculator555https://irsa.ipac.caltech.edu/applications/DUST. This service
provides position-based reddening estimates from Schlafly & Finkbeiner (2011).
Following Schlafly & Finkbeiner (2011), we use the Fitzpatrick (1999)
extinction curve for the Milky Way, with $R_{V}=A_{V}/E(B-V)=3.1$, where
$A_{V}$ is the extinction in magnitudes in the $V$ band, to calculate the
extinction at the central wavelength of each HST filter.

Finally, we perform a photometric correction on a small subsample of sources
in the catalog to account for elliptical apertures that were stretched or
overly elongated. This aperture stretching can occur for faint objects with
very close, bright neighboring sources. In these cases, while the fluxes
measured in the small Kron apertures may be affected, the aperture
corrections, which rely on the larger Kron aperture, are particularly
unreliable. Following Finkelstein et al. (2015b), we identified sources in the
catalog that needed this correction as those with a
$f_{\mathrm{AUTO}}/f_{\mathrm{APER}}$ flux ratio larger than the 95th
percentile of all sources in their $H_{160}$ magnitude bin. Here,
$f_{\mathrm{AUTO}}$ refers to the flux measured in the elliptical Kron
apertures, and $f_{\mathrm{APER}}$ is that measured in the circular $r=0\farcs
2$ apertures. We also required that the sources had a neighbor within 1″. We
then determined the median aperture corrections required to go from the
$r=0\farcs 2$ $f_{\mathrm{APER}}$ fluxes to the aperture-corrected
$f_{\mathrm{AUTO}}$ fluxes in bins of $H_{160}$ magnitude. We applied this
correction to the $f_{\mathrm{APER}}$ fluxes and uncertainties of the affected
sources. For these sources only, we use the corrected circular aperture fluxes
(rather than the Kron fluxes) in measuring source colors. This correction
affected only one of the sources in our sample of high-redshift candidates
(see Section 4).

###  2.5 Stars

We next construct a sample of stars in each of the datasets. This sample will
serve two purposes: to provide an estimate of the half-light radius
($r_{1/2}$) of unresolved sources in these images for comparison with our
selected candidates (Section 3.7), and to construct an average point spread
function (PSF) for use in simulating sources for our completeness simulations
(Section 6.2). We begin by selecting sources that are $\geq$5″ from any
neighboring objects in the catalog and that have measured half-light radii in
$H_{160}$ in the range $1<r_{1/2}<2.5$ pixels. We also require $H_{160}$
magnitudes in the range $16<m_{H160}<21$, where the cut at the bright end is
determined to avoid saturated sources and the cut at the faint end is to allow
for a clean separation between point sources and extended sources. Next, we
visually inspect the $H_{160}$ imaging of all star candidates, rejecting any
that lie along an image edge, are heavily affected by bad pixels or other
detector artifacts, or appear extended or as two very close sources that have
been extracted as a single source. Through this selection, we identify a
median of 4 moderately bright, isolated stars per field, though the number of
stars meeting these criteria in a given field ranges from 0 to 15. In total,
we identify 539 stars, 253 in BoRG[$z$8] fields, 112 in HIPPIES GO 12286
fields, and 174 in WISP fields.

Ideally, we would create a composite PSF for each filter in each HST pointing,
constructed from stars that are located across the full field. In this way, we
would account for any variations in the telescope focus during each set of
observations that could result in field-to-field variations in the PSF. By
including stars across the full field, we would also account for the detector
location dependence of the PSF. However, as these HST surveys were designed as
extragalactic programs, they attempted to avoid stars as much as possible.
Many of the HST pointings contain only a handful of isolated stars in the
desired magnitude range, and some fields contain none. As there are not enough
stars to construct a PSF for each individual HST pointing, we instead use a
single PSF for all fields of a given dataset. To validate this approach, we
explored the field dependence of the PSF by creating individual stacks of
stars for each field containing at least three point sources. We constructed a
radial profile of each combined PSF by computing the azimuthally-averaged flux
in circular apertures of increasing radii and used it to measure the full
width at half maximum (FWHM). We found that while the FWHM of the PSF varies
from field to field, it does so by an acceptably small amount. For example,
the standard deviation of $H_{160}$ FWHMs measured for individual fields is
$\leq$0$\farcs$04 for all three datasets. We therefore create a combined PSF
for each dataset/filter combination that we use for all fields in that
dataset.

For each dataset/filter, we create a combined PSF in the following way. The
individual image postage stamps of each identified star are resampled by a
factor of 10, registered to a common centroid using subpixel shifts, and then
resampled back to their original pixel scale. We normalize each image by the
flux measured in a circular aperture of radius $r=0\farcs 2$. Because the
fields are observed at a range of roll angles, we rotate each stamp by a
randomly-assigned angle in order to randomize the location of the diffraction
spikes. Finally, we median combine all image stamps. Combining multiple stars
more fully samples the undersampled HST/WFC3 PSF, as together they provide
images of point sources with a range of subpixel centroids. We perform this
stacking on each of the three datasets separately, to preserve the pixel
scales and drizzling patterns. The combined PSFs for each filter and dataset
are displayed in Figure 4, where we show the full stamps for all filters and
the central cores of the PSFs in the $H_{160}$ images as insets.

Figure 4:  Median-combined PSFs for each filter and dataset. The image stamps
are 10″$\times$10″, displayed on a zscale interval and linear stretch with no
image smoothing or pixel interpolation. We also show the central 3″$\times$3″
of the $H_{160}$ PSFs in the inset panels to zoom in on the core of the PSF.
The insets are displayed with a min-max interval and a log stretch. In each
panel, we list the measured FWHM of the PSF profile as well as the number of
stars that were combined in the stack. We display the $V_{606}$ PSF rather
than the $I_{600}$ PSF for the BoRG[$z$8] dataset, as $\sim$2/3 of the
BoRG[$z$8] fields are covered by $V_{606}$. The $I_{600}$ FHWM is included in
Table 3.

We calculate the radial profile and FWHM as described above. Each PSF panel in
Figure 4 is labeled with the measured FWHM and the number of stars that are
included in the stack. We also list the measured FWHMs in Table 3.

Table 3: Empirical PSFs and Point Source $r_{1/2}$

Filter | $N_{\mathrm{stars}}$ | Measured FWHM |  $H_{160}$ $r_{1/2}$ | $\sigma_{r1/2}$
---|---|---|---|---
|  | (arcsec) | (pixels) | (pixels)
BoRG[$z$8]
$V_{606}$ | 188 | 0$\farcs$179 |  $\cdots$ |  $\cdots$
$I_{600}$ | 77 | 0$\farcs$157 |  $\cdots$ |  $\cdots$
$Y_{098}$ | 253 | 0$\farcs$195 |  $\cdots$ |  $\cdots$
$J_{125}$ | 253 | 0$\farcs$181 |  $\cdots$ |  $\cdots$
$H_{160}$ | 253 | 0$\farcs$188 | 1.69 | 0.20
HIPPIES GO 12286
$I_{600}$ | 106 | 0$\farcs$161 |  $\cdots$ |  $\cdots$
$Y_{098}$ | 112 | 0$\farcs$181 |  $\cdots$ |  $\cdots$
$J_{125}$ | 112 | 0$\farcs$182 |  $\cdots$ |  $\cdots$
$H_{160}$ | 112 | 0$\farcs$178 | 1.39 | 0.14
WISP
$V_{606}$ | 173 | 0$\farcs$166 |  $\cdots$ |  $\cdots$
$I_{814}$ | 173 | 0$\farcs$177 |  $\cdots$ |  $\cdots$
$J_{110}$ | 174 | 0$\farcs$206 |  $\cdots$ |  $\cdots$
$H_{160}$ | 174 | 0$\farcs$181 | 1.37 | 0.26

Note. — Our measured FWHMs of the empirical, pixelated PSFs are larger than
the pre-pixelated values reported in the WFC3 Instrument Handbook (hst-
docs.stsci.edu/wfc3ihb).

Finally, we use this sample of visually-vetted stars to measure a typical
$r_{1/2}$ for stars in each dataset. While Source Extractor provides a
measurement of stellarity that can be useful in separating point and extended
sources (CLASS_STAR), we find that it is most reliable for bright stars,
$m_{H160}\lesssim 17$. The stars in this sample with $H_{160}$ magnitudes
fainter than $\sim$17 have CLASS_STAR values that range from 1.0 (as expected
for point sources) all the way down to $\sim$0.2 (closer to that expected for
extended sources). We therefore choose to use the $r_{1/2}$ as measured in
$H_{160}$ with Source Extractor to classify unresolved sources. In Table 3, we
provide the median and standard deviation of $r_{1/2}$ for the stars in each
dataset, values we use in Section 3.7 to idenfity potential stars in our
sample of high-redshift candidates.

##  3 Sample Selection

We select candidate $z\sim 9-10$ galaxies through a combination of criteria
related to detection significance, photometric redshift fitting, and visual
inspection. Whereas many have selected high-redshift galaxies through strict
color cuts (e.g., Oesch et al., 2013; Bouwens et al., 2015, 2016; Trenti et
al., 2011; Bradley et al., 2012; Schmidt et al., 2014a; Bernard et al., 2016),
we use a more holistic approach by including all available HST and Spitzer
imaging in our photometric redshift selection (as in, e.g., McLure et al.,
2010; Finkelstein et al., 2015b, 2022; Bouwens et al., 2019, 2021; Bowler et
al., 2020; Rojas-Ruiz et al., 2020). As the Lyman break is the main spectral
feature within the HST filter coverage for galaxies at $z\gtrsim 8$, the
photometric redshift fitting is akin to a Lyman break selection. However, this
approach allows for the selection of sources that might lie right on the edge
or even outside of a selection window in color-color space. Yet in
constructing a potentially more complete sample, we must take extra care with
contaminants. We therefore aim to be conservative and as objective as possible
in our process for source vetting, removing subjective steps entirely where
possible and quantifying our subjectivity where not. We describe the aspects
of our high-redshift galaxy selection and source vetting in the following
sections.

###  3.1 Photometric Redshifts

We use the Easy and Accurate $z_{\mathrm{phot}}$ from Yale (EAZY version
2015-05-08; Brammer et al., 2008) software to calculate photometric redshifts
for each source, real and simulated. The EAZY code fits all photometric
measurements to a synthetic spectral energy distribution (SED) template via a
chi-squared minimization process. In deriving a best-fitting template, it can
consider linear combinations of templates from a user-supplied set, and the
resulting fits are therefore less dependent on the choice of input templates.
We use the 12 tweak_fsps_QSF_12_v3 EAZY templates that are based on the
Flexible Stellar Population Synthesis (FSPS) models (Conroy et al., 2009;
Conroy & Gunn, 2010) and account for systematic differences between observed
galaxy colors and the models. Following Finkelstein et al. (2022), we also add
the spectrum of Q2343-BX418, a young ($<$100 Myr) galaxy at $z=2.3$ with high
equivalent width nebular emission lines (Erb et al., 2010). We include this
low mass ($M*\sim 10^{9}M_{\odot}$), low metallicity ($Z\sim 1/6Z_{\odot}$),
blue (UV continuum slope $\beta=-2.1$) template as galaxies with blue colors
are expected at high redshift (e.g., Bouwens et al., 2009; Finkelstein et al.,
2012b). We also add a version of this spectrum with all Ly$\alpha$ emission
removed to approximate attenuation by a neutral intergalactic medium (IGM)
while still providing a template with strong optical emission lines that can
affect Spitzer/IRAC colors (e.g., [O III] and H$\beta$ for galaxies at $z\sim
8-9$). We allow EAZY to construct best-fitting templates through linear
combination pulling from all 14 of these input templates.

In fitting templates over a grid of redshifts, EAZY compares the input source
fluxes and flux errors with the synthetic fluxes of the templates convolved
through each filter. We use the total fluxes computed in the small elliptical
Kron apertures as the input source fluxes (see Section 2.4, 2.4.2), with flux
uncertainties computed as described in Section 2.4.1. We also add a minimum
fractional error of 0.05 to the flux uncertainties in all filters (EAZY
parameter SYS_ERR) to allow for systematic uncertainty in our flux
measurements. Template fitting is performed in the redshift range $0.01\leq
z\leq 12$ in steps of $\Delta z=0.01$. At each redshift step, EAZY applies IGM
attenuation following Madau (1995) and includes the absorption by the
Ly$\alpha$ forest and damped Ly$\alpha$ systems as prescribed by Inoue et al.
(2014). As the population of galaxies at redshifts $8-10$ is currently not
well-understood, we assume a flat luminosity prior. This approach has the
benefit of not biasing us against potential true high-redshift sources, yet
also means we are treating sources at all redshifts the same, i.e., a galaxy
at $z\sim 2$ is equally likely as one at $z\sim 10$. We must therefore take
extra care in considering the contamination fraction of lower-redshift
galaxies in our sample, which we discuss in Section 5.

###  3.2 Selection Criteria

We begin by imposing an initial set of detection criteria on the full Source
Extractor catalog using the S/N as measured in circular apertures of radius
$r=0\farcs 2$. These small apertures capture the light at the central core of
the source positions as detected in the $H_{160}$ image and therefore maximize
the measured S/N. First, we require that sources be detected in the $H_{160}$
band with a S/N${}_{H}>5$. As we are aiming to select $z\gtrsim 9$ galaxies,
we expect the sources to have “dropped out” of the $V$, $I$, and $Y$ filters.
We therefore only consider sources with a S/N $<2$ in all available filters
blueward of the $J$ band (see Figure 2), permitting a S/N $>1.5$ in at most
one of these filters. This criterion allows for a detection up to S/N $=2$ in
one, but not both, of the bluer filters, such that a coincidental $<2\sigma$
noise fluctuation in an optical filter or a partial detection in, e.g.,
$Y_{105}$ for $8.5\lesssim z\lesssim 9$ will not disqualify a source.

We next select sources based on their $H_{160}$ magnitude, requiring
$22<H_{160}\leq 26.5$. The cut at the bright end helps remove contaminants
such as stars and lower-redshift galaxies, as $H_{160}\sim 22$ corresponds to
an absolute UV magnitude of $M_{UV}\sim-25$ at $z\sim 9$. We choose the
magnitude limit at the faint end in an attempt to remove spurious sources. As
our analysis involves relatively shallow HST imaging (as compared to the
Hubble Ultra Deep Field or the Hubble Frontier Fields, for example), we focus
on bright high-redshift candidates. The majority of the fields in these
datasets666We note that $\sim$30 of the fields we include in our analysis have
$5\sigma$ depths brighter than 26.5 as presented in Table 2. However, these
depths are averaged over the entire field and are not always indicative of the
depth in all areas of the field. have $5\sigma$ depths $\gtrsim$26.5, and so
the likelihood of detecting noise as spurious sources increases for magnitudes
fainter than 26.5. We use the total fluxes measured in the elliptical
apertures for these magnitude cuts.

As an additional check against including spurious sources in our selection, we
impose a minimum size criterion. Given the similarities in the depth and noise
properties of the parallel imaging we consider in this paper and that of
Rojas-Ruiz et al. (2020), we similarly expect a half-light radius criterion to
remove the majority of hot pixels and cosmic rays from our sample. Following
the example of Rojas-Ruiz et al. (2020), we used a random forest algorithm to
obtain a quantitative size cut for each dataset. In order to construct a
training set, we performed an initial visual inspection of 2500 randomly-
selected sources from each dataset that satisfied the detection significance
criteria in the $H_{160}$ and all optical filters. Based on our visual
inspection of the imaging in all available filters, we classified each source
as ‘real’ or ‘spurious’. These classifications became the target values or
class labels input to the random forest. To this sample, we added an
equivalent number of randomly-selected simulated sources that satisfied the
S/N criteria, each of which were classified as ‘real’. For each dataset, we
thus had a sample of 5000 sources that we fed to a Random Forest Classifier
using the Python Scikit-learn package (Pedregosa et al., 2011). Though we
included information such as the $H_{160}$ magnitude, source isophotal area,
elongation, and surface brightness, we found that the half-light radius was
the most discerning parameter. Through this analysis, we found that the
majority of sources classified as ‘spurious’ could be removed by imposing a
half-light radius cut at $r_{1/2}$$>$1.3 pixels for BoRG[$z$8] fields and
$>$1.1 pixels for the HIPPIES GO 12286 and WISP datasets. This difference in
$r_{1/2}$ values is expected due to the different pixel scales of each dataset
and corresponds to $r_{1/2}>0\farcs 1$ for the BoRG[$z$8] (0$\farcs$08/pixel)
and HIPPIES GO 12286 (0$\farcs$10/pixel) fields and $r_{1/2}>0\farcs 14$ for
the WISP (0$\farcs$13) fields. The $r_{1/2}$ values are measured by Source
Extractor in the $H_{160}$. We do not include a maximum effective radius
criterion in our selection, as it has been noted that such a criterion can
remove real sources along with lower-redshift interlopers (e.g., Holwerda et
al., 2020).

Finally, we use the redshift probability distribution functions (PDFs, which
we denote $p(z)$), that EAZY calculates for each source to select a sample of
candidates with preferred high-redshift solutions. Specifically, we require
that $>$70% of the integrated redshift probability be at $z>8$. This criterion
ensures that the majority of the redshift probability resides at high redshift
with less than 30% contained in lower-redshift solutions. We do not enforce
any additional criteria related to the photometric redshift fits at this
stage, as the PDFs derived from measurements in only four filters provide
minimal constraining power. For example, the redshift probabilty distribution
functions are flat for the majority of our candidates at $z\gtrsim 10$, where
the Ly$\alpha$ break has redshifted almost entirely out of the $J$ band and
the source is detected in only a single filter. While we do not select sources
based on the $\chi^{2}$ of the best-fitting template, we explore the
$\Delta\chi^{2}$ of the best-fitting templates at high and low redshift in
Section 5.

In summary, in order for sources to be considered $z\sim 9-10$ candidates, we
require that they:

  * •

are detected in $H_{160}$ with (S/N)${}_{H}>5$;

  * •

have $H_{160}$ magnitudes in the range $22\leq H_{160}\leq 26.5$;

  * •

are undetected at the $2\sigma$ level in $V$, $I$, and $Y$, allowing a
S/N$>$1.5 in at most one of these filters;

  * •

have a half light radius of $r_{1/2}$$>$1.3 (1.1) pixels for BoRG[$z$8]
(HIPPIES GO 12286 and WISP)

  * •

have at least 70% of their redshift probability distribution function at
$z>8$.

We find 193 sources that pass these selection criteria: 116 from BoRG[$z$8]
fields, 48 from HIPPIES GO 12286, and 29 from WISP fields. In the following
section we screen this sample for cases of detector persistence from
previously observed bright targets, which could masquerade as high-redshift
galaxies.

###  3.3 Persistence

Image persistence occurs when a bright source saturates the detector and
leaves a residual charge that appears as a ghost image in subsequent
exposures. This situation is particularly problematic in fields for which the
$H_{160}$-band is observed before the $J$-band, as persistence fades with time
and so can mimic a $J$-dropout. The BoRG observing strategy was designed to
mitigate image persistence by observing with the $Y$ band followed by
$J_{125}$ and $H_{160}$ in every orbit (Trenti et al., 2011). Persistence is
therefore expected to impact both the dropout and detection filters such that
the residual ghost images will not be selected as dropout candidates. However,
the other datasets we include in our analysis did not always perform their
observations in this order. Even some of the HIPPIES observations from GO
11702 that are included in the BoRG[$z$8] dataset were not scheduled such that
the $Y$-band images could shield the detection filters. The WISP observing
strategy is particularly prone to self persistence, which occurs when the
offending bright targets are observed in the same visit as the affected image.
As a slitless spectroscopic survey, WISP observed fields in direct
imaging–grism pairs, with the $J_{125}$ and $H_{160}$ images used for source
extraction and wavelength calibration of the G102 and G141 grism observations.
The zeroth orders of bright stars are then common causes of persistence in the
direct images that follow grism exposures.

We therefore performed a detailed persistence check for all candidates in the
following way. For each candidate, we searched the MAST archive for
observations taken in the 24 hours prior to each of the individual exposures
(FLTs) that went into the $H_{160}$ mosaics. If any of those observations had
a count $>$100,000 e- within a 10$\times$10 pixel box centered on the position
of a given candidate, the candidate was removed from our sample. We note that
the use of a 10 pixel box is likely overly conservative, and that the pixel
with the highest counts need not be the central pixel. However, considering a
box rather than a few central pixels allows for offsets in the drizzle
solution, as the world coordinate solution in the undrizzled image (that has
not been distortion corrected) could be off by a few pixels at the position of
the candidates. We identify 10 candidates as heavily affected by persistence
and remove them from our sample.

Figure 5: Two examples of WISP sources identified as self persistence and
removed from our sample. First, in panel (a) we show a 50$\times$50 pixel
stamp of a candidate identified as persistence due to the high detector counts
in FLTs observed before the $H_{160}$ observations. In panel (b), we show
10$\times$10 pixel stamps of the four affected FLTs that contributed to the
$H_{160}$ mosaic. The cause of the persistence was the zeroth order of a
bright star observed with the G141 grism just prior to the $H_{160}$
exposures. We show a zoom-in of this zeroth order in panel (c) and a larger
portion of the G141 image in panel (d) with the zeroth and first orders
indicated. As can be seen in panel (e), the position of self persistence
resulting from this specific observing strategy will be offset by
$\sim$180-195 pixels from the offending bright source. We have identified and
rejected two additional WISP sources in our sample based on their position
relative to a bright star, one of which is shown in panel (f). These
candidates were not identified as persistence based on counts from previous
observations at the candidate position on the detector, indicating that even
this thorough check can miss cases of persistence.

However, given the WISP observing strategy, we find that even this aggressive
approach is likely to have missed some cases of self persistence in WISP
images – where zeroth orders appear as $J_{110}$-dropout candidates. For each
source, the grism dispersion determines the location of all spectral orders in
the grism exposures with respect to the position of the source in the direct
image. We can therefore expect persistence from G141 zeroth orders to exist in
the $H_{160}$ image $\sim$180-195 pixels (23$\farcs$4-25$\farcs$4) to the left
of a bright star, where the exact offset is location dependent. We illustrate
this situation in Figure 5. First, in panel 5a, we show an example of a WISP
candidate that was identified as persistence based on the counts in
previously-observed FLTs as described above. This persistence is caused by a
zeroth order in an earlier G141 observation that is shown in panels 5c and 5d.
We show the $H_{160}$ image of this candidate in panel 5e, where it can be
seen that this candidate is at the expected distance to the left of a bright
star. In Figure 5f, we show another WISP candidate that was not identified by
our persistence check described above (the maximum counts in all FLTs was
$<$8,000 e-), yet lies at the expected position. We have removed 2 such WISP
sources from our sample based on their positions relative to a very bright
source. We note that we are only able to identify these additional cases of
persistence because of the specifics of the WISP observing strategy, and that
similar unidentified cases of persistence may still contaminate the rest of
our sample. Detecting candidates in a second imaging band is the only way to
confirm such single-band detections are real sources. We show all 12 sources
that are rejected as persistence in Figure 18d in the appendix, and indicate
those that were identified based on their distance from a bright star.

After removing the 12 cases of persistence, there are 181 remaining
candidates. We further explore these sources in Section 3.4 with a visual
inspection process to remove spurious or misclassified sources.

###  3.4 Visual Inspection

Thus far, we have selected candidate $z\sim 9-10$ galaxies using only
quantitative criteria. However, we find that the majority of sources in our
sample at this stage are spurious detections, where artifacts such as
diffraction spikes and partially-removed satellite trails are selected as
high-redshift candidates. Given the wavelength dependence of the PSF,
diffraction spikes will be largest in the $H_{160}$ images and can therefore
appear as $J$-band dropouts. Similarly, a satellite trail present in the
$H_{160}$ image can also cause false dropout candidates, and many of the
fields were observed with only two or three $H_{160}$ exposures, too few to
fully remove a satellite trail through outlier rejection during mosaic
creation. Spurious sources are therefore detected in the increased noise of
these satellite remnants. Hot pixels and cosmic rays that impacted the
detector in a compact area can also masquerade as high-redshift sources. The
drizzling process and kernel smoothing that is applied by Source Extractor can
both act to spread out the light into nearby pixels, giving cosmic rays a more
extended appearance like that of real sources. As these pure parallel datasets
are undithered, such features will not all be removed by image stacking. While
bad pixels on the IR detector should appear in the same location in all IR
filters, cosmic rays in the $H_{160}$ images can appear as $J$-band dropouts
in our selection. A careful visual inspection is necessary to identify and
remove such spurious candidates and other sample contamination.

However, visual inspection is an inherently subjective process, and the effect
that visually rejected candidates have on the sample incompleteness is
typically not accounted for with completeness calculations. We therefore use
simulated sources (see Section 6.2) to quantify any biases that may enter our
sample selection through the visual inspection process. We compile a set of
362 sources, consisting of all selected real candidates and an equal number of
simulated candidates selected in the same fields and via the same criteria. We
randomize this set of real and simulated candidates and carefully inspect each
one, accepting or rejecting them without knowing which is simulated. We
discuss using the simulated sources to quantify the effect of the visual
inspection on sample incompleteness in Section 6.2.1. For the rest of this
section we refer to the inspection and rejection of the real candidates. We
create a set of categories for rejection, ensuring that we use a standardized
classification scheme in all visual inspections. The classification categories
we consider are as follows.

First, we consider whether the source in question is a spurious detection. As
mentioned above, diffraction spikes and satellite remnants are common examples
of spurious detections. Additionally, we reject sources located immediately
along image edges where the noise levels are significantly higher, because
these sources at best have unreliable photometry and at worst are only
partially extracted. We also expect spurious sources may be identified on the
edges of bright sources, where the Source Extractor deblending can be too
aggressive. In these oversplit regions, the light from a bright source gets
separated and attributed to multiple sources. The PSF is smaller in bluer
filters, and so any oversplit regions are closer to the bright object core.
Oversplit regions in $J$ and $H$ can therefore appear as high-redshift
galaxies that have dropped out of the bluer filters. We consider this
possibility during our inspection, but do not identify any sources in this
category. As part of our visual inspection, we also inspect the RMS maps at
the position of each candidate and reject sources with clusters of bad pixels
identified within the $r=0\farcs 2$ circular apertures used to determine
detection significance. Our final category of spurious detections includes hot
pixels, cosmic rays, and other artifacts that are not identified as bad pixels
in the RMS maps. We identify these features visually as sources appearing as
single pixels, disjointed clusters of pixels that may have been smoothed
together during source detection, or otherwise strange morphologies.

We then consider two additional categories for candidate rejection, cases in
which a source is identified as real yet was incorrectly selected due to
problems related to photometric measurement or redshift fitting. As described
in Section 2.4.2, the Source Extractor Kron apertures of small sources in
close proximity to a large, bright neighbor can be stretched by the light from
the neighbor. We have attemped to identify sources with affected apertures and
correct the corresponding fluxes (see Section 2.4.2), and one of the 181
selected candidates required this correction (Par0456-2203_473, see Section
4.1.4). Here, we use our visual inspection process to check that the
elliptical apertures (or, in the case of the one corrected source, the
circular aperture that replaced the elliptical aperture) closely match the
source morphology for all candidates. We also inspect all available $V$, $I$,
and $Y$ band images at each source position to ensure that there is no
significant flux that was missed in the catalog measurements.

In summary, we consider the following categories of spurious sources or
unreliable measurements in our visual inspection (and indicate in parentheses
the number of real sources rejected in each category):

  1. 1.

Is the candidate a real source or a spurious detection (166 rejected, see
Appendix B)? I.e., is the source:

     * 1a.

a diffraction spike or satellite remnant (48 rejected, see Figure 18a);

     * 1b.

along an image edge or detected in an area of increased noise near an image
edge (12, Figure 18b);

     * 1c.

a bad pixel in the weight map (25, Figure 18c);

     * 1d.

the oversplit region of a bright neighboring source (0); or

     * 1e.

a hot pixel, cosmic ray, or source with a strange morphology (81, Figure 19)?

  2. 2.

Is the aperture drawn correctly, or was the Kron ellipse stretched by a
neighboring source and not subsequently corrected? (0)

  3. 3.

Is there significant optical flux that was not measured correctly in the
catalog? (0)

As can be seen, all 166 candidates we rejected through visual inspection were
identified as some type of spurious detection or false source. While 85 of
these rejections are objectively motivated (diffraction spikes, satellite
trails, image edges, and bad pixel clusters), 81 are based on subjective
opinion during the inspection (category 1e). In Section 6.2.1, we use the
classifications of simulated sources performed as part of the “anonymous”
visual inspection to explore these 81 classifications.

Following this visual inspection, there are 15 candidates that have satisified
all of our selection criteria. We further explore these sources in Sections
3.5, 3.7, and 5.

###  3.5 IRAC

With 15 candidates remaining after rejecting sources from our visual
inspection and persistence checks, we now turn to incorporating Spitzer/IRAC
(Fazio et al., 2004) imaging to our analysis. Imaging at these redder
wavelengths is a crucial tool for disentangling the spectral energy
distributions of true high-redshift galaxies and sources that are common
contaminants in high-redshift samples (Finkelstein et al., 2022, e.g.,). The
$J-H$ colors of both passive and dusty star-forming galaxies at $z\sim 2-3$
can be indistinguishable from those at $z>8$, and these contaminants will
likely by undetected in bluer filters given the shallow imaging of these
parallel datasets. Additionally, as the majority of galaxies at $z>8$ are
compact or unresolved in the HST NIR imaging, M, L, and T dwarf stars can
similarly contaminate HST-selected high-redshift samples. However, the
spectral energy distributions of all three types of contaminants are expected
to diverge from those of high-redshift galaxes at $\lambda\gtrsim 2$µm, and
even shallow IRAC imaging can help distinguish between them.

Table 4: IRAC Coverage  Field | Program | PI | Ch1 AORKEY | Ch2 AORKEY
---|---|---|---|---
Par335 | 10041 | Colbert | 51773440 |  $\cdots$
Par0440-5244 | 14253 | Stefanon | 69070080 | 69070336
Par0456-2203 | 12058 | Bouwens | 59360512, 59360768 |  $\cdots$
Par0456-2203 | 14253 | Stefanon |  $\cdots$ | 69071360
Par0713+7405 | 11121 | Finkelstein | 53188864 |  $\cdots$
Par0713+7405 | 12058 | Bouwens | 58126848, 58127104 |  $\cdots$
Par0713+7405 | 14253 | Stefanon | 69082368 | 69075200, 69082624
Par0756+3043 | 12058 | Bouwens | 58126336, 58126592 |  $\cdots$
Par0756+3043 | 14253 | Stefanon |  $\cdots$ | 69078016
Par0843+4114 | 14253 | Stefanon | 69082368 | 69082624
Par0926+4000 | 14304 | Stefanon | 69982720 |  $\cdots$
Par0926+4536 | 11121 | Finkelstein | 53189120, 53189376 |  $\cdots$
Par0926+4536 | 12058 | Bouwens | 58125568, 58125824 |  $\cdots$
Par0926+4536 | 14253 | Stefanon |  $\cdots$ | 69091328
Par0956-0450 | 11121 | Finkelstein | 53189632 |  $\cdots$
Par0956-0450 | 14253 | Stefanon |  $\cdots$ | 69096704
Par1033+5051 | 12058 | Bouwens | 58125056, 58125312 |  $\cdots$
Par1033+5051 | 14253 | Stefanon | 69105152, 69105408 |  $\cdots$
Par1301+0000 | 80134 | Colbert | 42568960 |  $\cdots$
Par1301+0000 | 14253 | Stefanon | 69136896 | 69137152
Par2346-0021 | 90045 | Richards | 46957568, 46961152 | 46957568, 46961152

We searched the Spitzer Heritage Archive hosted by
IRSA777sha.ipac.caltech.edu/applications/Spitzer/SHA, Spitzer Level 2 data:
www.ipac.caltech.edu/doi/irsa/10.26131/IRSA413 for all available IRAC imaging
at the positions of each candidate and downloaded the coresponding Level 2
(“Post Basic Calibrated Data”; PBCD) mosaic images. These mosaics are
processed by the IRAC pipeline version S19.2888except for Program 12058
coverage of Par0926+4536 and Par1033+5051, which used version S.19.1. See
irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook for
information about the IRAC pipeline. and are on a 0$\farcs$6 pixel scale.
Spitzer/IRAC 3.6µm imaging exists for 11 of the HST parallel fields in which
we identify candidates (nine of which also have 4.5µm imaging), amounting to
IRAC coverage of all but two of the 15 candidates. The imaging was obtained as
part of seven unique programs, many of which were designed with the goal of
following up high-redshift candidates previously identified in these fields.
Table 4 lists the program information for each IRAC dataset that we include in
our analysis.

We measure IRAC photometry using a method similar to that described in Rojas-
Ruiz et al. (2020). For each candidate, we create background-subtracted IRAC
stamps that are 30$\farcs$6$\times$30$\farcs$6 (51$\times$51 pixels). We
consider each IRAC observation separately, such that a source observed at
multiple position angles or through multiple programs may have multiple stamps
per channel. This approach results in 36 stamps (24 of which are at 3.6µm) for
the 13 sources with IRAC coverage. In our photometric analysis, we treat each
stamp as an independent measurement. We model all sources in each stamp with
the Galfit999users.obs.carnegiescience.edu/peng/work/galfit/galfit.html image-
fitting software (v3.0 Peng et al., 2010) as demonstrated by Finkelstein et
al. (2015a). Given the lower resolution of IRAC images, this modeling approach
is necessary to deblend the light from the candidate and any neighboring
sources. We use source positions and magnitudes from the $H_{160}$ catalog as
inputs to Galfit, including everything in the catalog down to $H_{160}=25$ and
making the 3.6µm and 4.5µm magnitude initial guesses one magnitude brighter
than the $H_{160}$ value. We constrain source positions in the Galfit model to
be within $\pm$1.5 pixels (0.9″) of the input values and model magnitudes to
be brighter than 40. We can safely assume that sources that hit this lower
magnitude constraint are undetected in the IRAC imaging, and so we iteratively
remove them from the model. Extended sources, defined as those with a semi-
major axis in the $H_{160}$ Source Extractor catalog that is $>$2$\times$
larger then the FWHM of the IRAC point response function, are modelled as
Sérsic profiles. However, the majority of sources in the IRAC stamps – all
high-redshift candidates and almost all of their neightbors – are modelled as
point sources.

Figure 6: The median-combined PSFs created from IRAC 3.6µm point sources (left
two columns) and 4.5µm point sources (right two columns). The separate PSFs
for Par2346-0021 are labeled ‘Par2346’. The top row displays the full stamps
($30\arcsec\times 30\arcsec$) on a zscale interval and linear stretch, and the
number of individual point sources that were combined to create each PSF are
indicated by $N$. The bottom row zooms in on the central 15″ and uses a min-
max interval with a square root stretch to highlight the PSF core.  Figure 7:
Spitzer/IRAC imaging of the 13 candidates. For each candidate we show a
30$\farcs$6 IRAC image centered at the candidate’s position (left), the
30$\farcs$6 residual map with all source flux removed except that of the
candidate (middle), and a 10$\farcs$2 zoom-in of the residual map to highlight
the structure close to the candidate position (right). In each set of stamps,
the top and bottom rows display the 3.6µm and 4.5µm (where available) images,
respectively. The black squares indicate the size of the zoom-in stamp, and
the circles have radius 2.5 pixels, the size used for aperture photometry as
described in the text. For fields observed multiple times, we show a single
observation as an example. Additional images and residual maps are similar to
those shown here. The image stamps are displayed with no image smoothing and
on zscale intervals calculated individually for each stamp and zoom. The
deblending at the target positions is mostly clean with the exception of
Par0456-2203_473, which is discussed in Section 4.1.4.

For the point source models, we create median PSFs for each IRAC channel in
the same manner as described in Section 2.5. Specifically, we run Source
Extractor on the IRAC images and identify point sources as bright
($14<M_{\mathrm{IRAC}}<18$), isolated (no neighbors within 15″) sources with a
measured half light radius between
$0.9\arcsec<r_{1/2,\mathrm{IRAC}}<1.35\arcsec$. We visually inspect all
sources to remove any with bad pixels or undetected neighbors in the Source
Extractor catalog, resample all image stamps by a factor of 10 to allow for
sub-pixel centroiding, and median combine the resulting stack. We do not
randomly rotate each stamp, as all the IRAC mosaics are aligned by detector
coordinates. In total, we combine 243 (3.6µm) and 88 (4.5µm) sources to create
the PSFs, which are shown in Figure 6. In Figure 7, we show the results of
running Galfit on the 13 candidates with IRAC coverage.

We use a two-fold approach to measuring the IRAC photometry of the sources in
our sample. First, we adopt the Galfit model magnitudes for all IRAC
observations in which Galfit successfully modeled the target source, as was
the case for 23 out of the 36 separate stamps. For the remaining 13, Galfit
did not measure a significant flux at the source position, i.e., the Galfit
model magnitude hit the constraint at 40 mag and the source was removed from
the model. In these cases, we measure the photometry of the target in the
Galfit residual maps with all neighboring source flux removed. We use circular
apertures of radius $r=2.5$ pixels, corresponding to a diameter of 3″
($\sim$1.5$\times$ the FWHM of the 3.6µm warm mission point response
function101010https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/iracinstrumenthandbook/5).
The black circles in the zoomed-in residual maps in Figure 7 indicate the
aperture sizes used to measure the candidate fluxes. We apply aperture
corrections of 1.69 (3.6µm) and 1.70 (4.5µm), obtained by measuring the flux
of our custom PSFs in circular apertures of increasing radii.

We estimate the flux uncertainty in these measurements in the follow way. We
place $r=2.5$ pixel circular apertures randomly across the full, background-
subtracted IRAC images, using Source Extractor segmentation maps to avoid
source flux. With Photutils, we measure the flux in these apertures, fit a
Gaussian to the flux distribution, and take the standard deviation as the flux
uncertainty in the apertures. We then treat this $1\sigma$ uncertainty as a
noise floor, adopting it as the minimum uncertainty for both the Galfit model
fluxes and those measured through aperture photometry. Indeed, the majority of
the Galfit uncertainties are larger than the background-measured standard
deviation, indicating that the Galfit photometry is also accounting for
uncertainties related to the source flux deblending.

Finally, we rerun EAZY for each high-redshift candidate, incorporating all
available HST and Spitzer/IRAC photometry. As mentioned previously, we fit
each IRAC observation with Galfit separately, leading to multiple photometric
measurements for some candidates. We adopt a weighted mean of the available
measurements for each IRAC channel in our EAZY runs. We note that the EAZY
results – including the best-fitting SED templates and the redshift
probabiltiy distributions – are nearly identical if we instead treat each IRAC
flux as an independent measurement, i.e., with multiple instances of the IRAC
$3.6\mu$m and 4.5µm filters in the EAZY setup files. The best-fitting
redshifts from each method (weighted mean versus multiple photometric
measurements) differ by at most $\Delta z<0.05$.

With the exception of one candidate (Par2346-0021_164, discussed in Section
3.5.1), the HST and Spitzer/IRAC photometry continue to prefer a high-redshift
solution. For many of these candidates, the addition of the IRAC photometry
reduces the fraction of the redshift PDF that lies at $z<6$, decreasing this
fraction from $\sim$0.1-0.2 in most cases to $<$0.01. The most extreme
improvements were for candidates Par0926+4000_369
($p(z$$<$$6)_{\mathrm{HST}}=0.295$ to $p(z$$<$$6)_{\mathrm{HST+IRAC}}<0.001$)
and Par0843+4114_120 ($p(z$$<$$6)_{\mathrm{HST}}=0.236$ to
$p(z$$<$$6)_{\mathrm{HST+IRAC}}=0.009$). This rejection of the lower-redshift
solutions typically occurs when the upper limits measured in one or both of
the IRAC channels fall below the flux density expected for lower-redshift, red
and dusty galaxies. For some candidates, the EAZY run with IRAC photometry
even helped narrow the redshift probability peak at $z>8$, thereby tightening
the photometric redshift constraints. The 95% interval of the redshift PDF at
$z>8$ decreases by $z\sim 0.36$ on average and as much as 0.9 for candidate
Par0926+4000_369. In these cases the IRAC measurements, as either detections
or upper limits, help constrain the expected location of the 4000Å break. We
show the best-fitting SED templates and photometric redshift PDFs
incorporating HST and IRAC photometry for our high-redshift candidates in
Figures 10 and 11, and discuss the IRAC results of individual candidates in
Section 4.1.

####  3.5.1 Field Par2346-0021

While almost all of the IRAC observations listed in Table 4 covered a single
pointing with relatively small dithers, those covering HIPPIES field
Par2346-0021 (PI: Richards; PID 90045) are a little different. This program
aimed to cover a large enough volume to study the clustering and luminosity
function of quasars at $z>3$, and so it created wide-area maps along the Sloan
Digital Sky Survey (SDSS) ‘Stripe 82’. The PBCD maps at the position of our
high-redshift candidate cover $\sim$$2.2^{\circ}\times 0.65^{\circ}$ and are
significantly shallower than the IRAC imaging available for the rest of our
sample. The 1$\sigma$ flux uncertainties in these large maps, measured in
$r=2.5$ pixel circular apertures as described above, are $\sim$1545 nJy and
$\sim$1407 nJy at 3.6µm and 4.5µm, respectively, while the median
uncertainties for the rest of the fields are 302 nJy and 737 nJy. This
difference corresponds to limiting magnitudes that are $\sim$1.8 (3.6µm) and
0.7 (4.5µm) magnitudes brighter in this field than the others we consider in
Section 3.5. The much larger imaging area also provides far more point sources
for the creation of a custom PSF.

We therefore measure IRAC photometry for the candidate in Par2346-0021
separately, following the same steps described in Section 3.5 but with field-
specific custom PSFs. We create the PSFs for this field using only sources
identified in the Par2346-0021 IRAC maps, though we use a slightly more
conservative magnitude cut when identifying point sources because there are so
many to work with: $14<M_{\mathrm{IRAC}}<17.5$. The PSFs we create for
Par2346-0021 are comprised of 1561 (3.6µm) and 1259 (4.5µm) sources and are
also displayed in Figure 6 in columns titled ‘Par2346’. The aperture
corrections we measure using the 3.6µm PSF (1.67) is similar to that measured
for the other IRAC observations, and the 4.5µm aperture correction (1.98) is
slightly larger. We then proceed as before with Galfit and measuring source
photometry and flux uncertainties.

The Galfit model brightnesses calculated for Par2346-0021_164, when combined
with the HST photometry, reduce our confidence in this candidate. As shown in
Figure 8, this source has a very red $H-[3.6]$ color that is consistent with a
dusty, red SED. The photometric redshift PDF obtained with the HST+IRAC
photometry (blue, solid line) has an increased peak around $z\sim 3$ when
compared with that from the HST photometry alone (purple, dashed). Due to this
added probability at lower redshifts, only 62% of the integrated PDF is at
$z>8$, rather than the 70% threshold we require as part of our selection
criteria.

While we may opt to keep an HST-selected source in our sample if we find that
we cannot trust the IRAC photometry, that is not the case for this candidate.
It lies outside of the existing $I_{600}$ imaging in this field, and so is
only covered by three HST filters, and the HST-only photometric redshift
solution is not well-constrained (the 68% redshift interval ranges from
$z=3.07$ to 10.06). Additionally, the candidate does not appear to have any
close neighbors in the $H_{160}$ image, and so we consider the IRAC photometry
reliable. Although the IRAC imaging is shallow, the Galfit model magnitudes
for the candidate are significant. We therefore remove Par2346-0021_164 from
our sample as a possible $z<7$ contaminant, but note that it could be a viable
$z\sim 9$ candidate.

Figure 8: The image stamps, photometry, best-fitting SEDs and redshift PDFs
for two candidates we have removed from our sample. For each candidate we show
3″ HST stamps and, where available, 10″ Galfit residual maps for the IRAC
imaging from each individual AOR with all modeled flux removed except that
attributed to the candidate. The purple $r=0\farcs 5$ circles identify the
position of the candidate, and the blue 3″ squares indicate the size of the
HST stamps on the larger IRAC images. Beneath the image stamps we show the
SEDs in the left panel and the photometric redshift PDFs in the right, with
the fit to HST+IRAC photometry (solid blue curve) and HST-only photometry
(dashed purple curve). Top: We remove Par2346-0021_164 due to an increase in
its redshift PDF at $z\sim 2-4$ after including IRAC photometry. Bottom: We
remove Par2132+1004_509 from our sample primarily due to its low
$\Delta\chi^{2}=2.53$, but also because it is only observed in 3 HST bands and
is not covered by Spitzer/IRAC.

###  3.6 Photometric Redshift Goodness-of-Fit

While we did not include a minimum $\chi^{2}$ threshold on the EAZY templates
as part of our selection criteria, we now explore the goodness-of-fit of these
photometric redshifts, including IRAC where available. Specifically, we
consider the difference between the $\chi^{2}$ of the best-fitting template
and the minimum $\chi^{2}$ at $z<6$. This $\Delta\chi^{2}$ provides an
estimate of the quality of the photometric redshift fit at high-redshift
compared to the best-possible fit at lower redshift. In all cases, the
$\chi^{2}$ of the best-fitting template is lower than that at $z<6$, with a
median $\Delta\chi^{2}=\chi^{2}_{z<6}-\chi^{2}_{\mathrm{best}}=16.8$. However,
the $\Delta\chi^{2}$ for candidate Par2132+1004_509 is 2.53, lower than the
value corresponding to the 95% confidence interval ($\chi^{2}>4$). This
candidate also has a broad PDF (the redshift interval containing 68% of the
PDF ranges from $z=4.01$ to 11.51), is not covered by IRAC imaging, and falls
on the UVIS chip gap in the $I_{600}$ imaging (see Figure 8). We therefore
consider this source as an unreliable high-redshift candidate and remove it
from our sample.

###  3.7 Comparison with Stellar Colors

Finally, we explore whether any of the $z\gtrsim 9$ candidates are likely to
in fact be stars. Low-mass stars and brown dwarfs can have NIR colors that are
very similar to those of high-redshift galaxies, and so are common sources of
contamination in Lyman break samples. The case for potential contamination is
made worse by the fact that all but one of the candidates in our sample are
unresolved in the HST imaging. Eight of the candidates have half-light radii
($r_{1/2}$) as measured in $H_{160}$ that are smaller than the median
$r_{1/2}$ of the stars in their respective datasets (see Section 2.5). The
$r_{1/2}$ of another four candidates are within $1\sigma$ of the median
values. The one exception is Par0456-2203_473, which has a larger $r_{1/2}$
and is discussed further in Section 6.1. The very compact sizes measured for
the majority of our sample are not necessarily concerning. These sources have
fairly low S/N in undithered data, and so are not expected to show significant
extended morphologies. However, it does make a comparison with stellar colors
crucial.

We therefore compare the colors of the 13 remaining candidates with those of
M, L, and T dwarfs. The stellar colors are calculated from spectra in the
range $\sim$0.63$-$2.5µm taken with the medium-resolution ($R\sim 2000$)
spectrograph SpeX (Rayner et al., 2003) at the NASA Infrared Telescope
Facility (IRTF). We downloaded from the IRTF Spectral
Library111111irtfweb.ifa.hawaii.edu/~spex/IRTF_Spectral_Library (Burgasser,
2014) 582 observations of 449 unique stars (132 M dwarfs, 210 L dwarfs, and
107 T dwarfs), covering the full range of temperature subclasses and including
$J$, $H$, and $Ks$ photometry from the Two Micron All Sky Survey (2MASS).

Figure 9: The SpeX-derived goodness-of-fit $\chi^{2}$ values as a function of
spectral type (orange, top axis) compared with that determined by EAZY as a
function of redshift (blue, bottom axis) for high-redshift candidate
Par0713$+$7405 ID 95. The $\chi^{2}$ calculated for each SpeX spectrum is
identified by a separate point, indicating the variation measured for a given
spectral type. The minimum galactic and stellar $\chi^{2}$ values are
identified by horizontal dashed and dashed-dotted lines, respectively. The
difference in $\chi^{2}$ values for this source is
$\Delta\chi^{2}=\chi^{2}_{\mathrm{stellar}}-\chi^{2}_{\mathrm{galactic}}=11.11$,
which is the smallest $\Delta\chi^{2}$ of all candidates in our sample. All
candidates are better fit by galactic templates than by MLT dwarf spectral
templates.

We calculate broadband HST fluxes by integrating each SpeX spectrum through
the HST filter profiles shown in Figure 2. The filters are interpolated onto
the observed wavelength array of each spectrum such that only the portion of
the filter profile that overlaps the spectrum is used in calculating broadband
fluxes. For each SpeX spectrum, we remove any fluxes that result from a $<$25%
wavelength overlap between the filter and the spectrum. Following Finkelstein
et al. (2022), we then assign IRAC magnitudes to each stellar spectrum using
the median $Ks-[3.6]$ and $[3.6]-[4.5]$ presented by Patten et al. (2006). We
match the SpeX stars to the colors in Tables 1 and 3 from Patten et al. (2006)
by spectral type, using $\Delta\mathrm{type}<1$ where possible and
$\Delta\mathrm{type}<3$ for a few sources with no closer spectral matches. We
convert all photometry – 2MASS $J$, $H$ and $Ks$ from the SpeX Library, and
3.6µm and 4.5µm obtained by matching to Patten et al. (2006) – from Vega to AB
magnitudes using $m_{\mathrm{Vega}}+0.91$, 1.39, 1.85, 2.79 and 3.26 for 2MASS
$J$, $H$, $Ks$, and Spitzer/IRAC 3.6µm and 4.5µm, respectively.

Finally, with a full suite of photometry for each SpeX spectrum, we calculate
a $\chi^{2}$ goodness of fit for the photometry of each high-redshift
candidate. We then compare this SpeX-derived $\chi^{2}$ with that from the
best-fitting galactic template determined by EAZY and find that
$\chi^{2}_{\mathrm{EAZY}}<\chi^{2}_{\mathrm{SpeX}}$ for all candidates. As an
example, in Figure 9, we show the $\chi^{2}$ measured for each SpeX specturm
given the photometry of Par0713$+$7405_95. Of all candidates in our sample,
this source has both the lowest minimum SpeX-derived $\chi^{2}$ (=13.07) and
the smallest $\Delta\chi^{2}(\mathrm{SpeX}-EAZY)=11.1$.

Therefore, although almost all of the 13 candidates are unresolved, none of
their HST/WFC3 and Spitzer/IRAC (when available) colors are better fit by
these 500$+$ M, L, and T dwarf spectra. We conclude that all 13 candidates are
better fit by galactic spectral templates and do not remove any from the
sample due to stellar contamination. Future observations with JWST/NIRCam
imaging, for example, will vastly improve on the stellar contamination
analysis currently possible in high-redshift Lyman break samples by providing
both higher-resolution imaging in the NIR and imaging in 5$+$ filters redward
of the expected position of the Lyman break.

Table 5: Sample of $z\sim 9-10$ Candidates

Par | ID | Dataset | RA | Dec | $m_{160}$ | $r_{1/2}$ | $z_{\mathrm{phot}}$ | $p(z>8)$ | 68% | $\Delta\chi^{2}$ | $M_{\mathrm{UV}}$
---|---|---|---|---|---|---|---|---|---|---|---
Par1033+5051 | 116 | BoRG[$z$8] | 10:32:44.7 |  $+$50:50:30.1 | 25.99 | 0$\farcs$11 | 8.28 | 81.29% | 7.86$-$8 | 21.25 | $-21.20^{+0.15}_{-0.14}$ |  |
Par0440-5244 | 497 | BoRG[$z$8] | 04:39:47.0 |  $-$52:43:55.2 | 25.47 | 0$\farcs$15 | 8.53 | 89.06% | 7.98$-$9 | 23.59 | $-21.76^{+0.17}_{-0.15}$ |  |
Par0713+7405 | 95ccThese candidates were found to experience intermediate-to-strong lensing, and so the $M_{\mathrm{UV}}$ values have been corrected for magnification as described in Section 6.1. | HIPPIES 12286 | 07:13:13.2 |  $+$74:05:04.5 | 25.74 | 0$\farcs$14 | 9.23 | 98.16% | 8.79$-$9 | 7.33 | $-21.32^{+0.18}_{-0.19}$ |  |
Par0456-2203 | 473ccThese candidates were found to experience intermediate-to-strong lensing, and so the $M_{\mathrm{UV}}$ values have been corrected for magnification as described in Section 6.1. | BoRG[$z$8] | 04:55:51.6 |  $-$22:01:36.1 | 24.52 | 0$\farcs$33 | 9.39 | 93.45% | 8.81$-$9 | 4.70 | $-19.87^{+0.56}_{-0.36}$ |  |
Par0259+0032aaThis field does not have Spitzer/IRAC coverage. | 194ccThese candidates were found to experience intermediate-to-strong lensing, and so the $M_{\mathrm{UV}}$ values have been corrected for magnification as described in Section 6.1. | HIPPIES 12286 | 02:59:37.6 |  $+$00:33:01.8 | 24.79 | 0$\farcs$13 | 9.36bbThis redshift is based on an upper limit estimated for a contaminated $J_{125}$ flux. See text for details. | 98.71% | 8.77$-$9 | 8.81 | $-21.13^{+1.06}_{-0.59}$ |  |
Par0926+4000 | 369 | BoRG[$z$8] | 09:25:33.6 |  $+$40:00:44.9 | 24.68 | 0$\farcs$13 | 10.49 |  $>$99.99% | 10.10$-$11 | 28.66 | $-22.88^{+0.11}_{-0.12}$ |  |
Par335 | 251 | WISP | 15:47:43.4 |  $+$20:58:23.3 | 24.74 | 0$\farcs$19 | 10.53 |  $>$99.99% | 10.24$-$11 | 20.41 | $-22.85^{+0.12}_{-0.12}$ |  |
Par0956-0450 | 684 | HIPPIES 12286 | 09:56:50.7 |  $-$04:49:38.8 | 24.88 | 0$\farcs$12 | 10.54 |  $>$99.99% | 10.49$-$11 | 69.05 | $-22.74^{+0.08}_{-0.09}$ |  |
Par0843+4114 | 120 | HIPPIES 12286 | 08:43:25.9 |  $+$41:14:00.3 | 26.39 | 0$\farcs$13 | 10.98 | 99.07% | 10.30$-$11 | 7.79 | $-21.23^{+0.24}_{-0.20}$ |  |
Par1301+0000 | 37 | BoRG[$z$8] | 13:01:13.6 |  $-$00:01:26.6 | 25.42 | 0$\farcs$13 | 11.02 | 99.99% | 10.46$-$11 | 16.49 | $-22.20^{+0.17}_{-0.15}$ |  |
Par0756+3043 | 92 | BoRG[$z$8] | 07:55:54.4 |  $+$30:42:22.1 | 26.03 | 0$\farcs$14 | 11.03 | 99.83% | 10.39$-$11 | 10.83 | $-21.59^{+0.22}_{-0.20}$ |  |
Par0440-5244 | 593 | BoRG[$z$8] | 04:39:55.2 |  $-$52:43:49.6 | 26.46 | 0$\farcs$11 | 11.07 |  $>$99.99% | 10.46$-$11 | 17.99 | $-21.16^{+0.15}_{-0.14}$ |  |
Par0926+4536 | 155 | HIPPIES 12286 | 09:26:26.2 |  $+$45:37:26.2 | 26.37 | 0$\farcs$13 | 11.15 | 99.99% | 10.59$-$11 | 16.78 | $-21.26^{+0.23}_{-0.20}$ |  |

Note. — Columns are: (1-2) field and candidate IDs; (3) Dataset containing
field; (4-5) candidate right ascension and declination; (6) $H_{160}$ apparent
magnitude; (7) half-light radius as measured in $H_{160}$; (8) best-fitting
redshift from EAZY, selected as the peak of the redshift PDF; (9) percentage
of the integrated PDF that is at $z>8$; (10) 68% interval of the redshift PDF;
(11) difference in the $\chi^{2}$ of the best-fitting EAZY template and the
minimum $\chi^{2}$ at $z<7$:
$\Delta\chi^{2}=\chi^{2}_{z<7}-\chi^{2}_{\mathrm{best}}$; (12) median absolute
magnitude and 68% interval of magnitude distributions as described in Section
6.1.

##  4 Results

###  4.1 High-redshift Candidates

After applying all selection criteria, visual inspection, persistence
screening and stellar color comparison, we have a sample of 13 high-redshift
candidate galaxies spanning the range $8.3\lesssim z\lesssim 11$. In Table 5,
we list the candidates, their positions, $H_{160}$ magnitudes, half-light
radii measured in $H_{160}$, and their photometric redshifts. The photometric
redshift fits all include IRAC where available. We also provide the percentage
of the integrated photometric redshift PDF that lies at $z>8$ (recall that we
require $p(z>8)>70\%$ for sample selection) and the redshift interval
containing 68% of the redshift PDF in columns (9) and (10), respectively. The
lower bound of this interval is above $z=8$ for all but two of the candidates,
and these both have $p(z>8)>80\%$. Column (9) also demonstrates that while the
redshift PDFs are relatively well-constrained to be at high redshifts, they
are still quite broad with 68% $p(z)$ intervals ranging from $\sim$0.8 to
$\sim$1.4. Column (11) lists the difference between the $\chi^{2}$ of the
best-fitting template and the minimum $\chi^{2}$ at $z<7$ for the same EAZY
run (see discussion in Section 3.5).

Figure 10: The image stamps, photometry, best-fitting SEDs and redshift PDFs
for the first 6 candidates in our sample, sorted by photometric redshift. For
each source, we show all available imaging including 3″ HST stamps and 10″
Galfit residual maps for the IRAC imaging from individual AORs with all
modeled flux removed except that attributed to the candidate. The purple
$r=0\farcs 5$ circles identify the position of the candidate in the HST
images, and the blue 3″ squares indicate the size of the HST stamps in the
larger IRAC images. Beneath the image stamps we show the SEDs in the left
panel and the photometric redshift PDFs in the right. We show the measured
photometry as filled circles and the predicted fluxes for each template as
open squares. The blue solid curve and blue shaded regions indicate the best-
fitting template and PDF including all HST+IRAC photometry. The dashed purple
curves indicate the best-fitting template/PDF considering only HST photometry.
The inclusion of the IRAC photometry almost always tightens the width of the
PDF at high redshifts while reducing the portion of the PDF present at lower
redshifts. The red dot-dashed curve shows the best-fitting template/PDF at low
redshift, when EAZY is restricted to $z<6$. These candidates are all better
fit with a high-redshift galactic spectral template than that of a lower-
redshift galaxy. We use a different color scheme with Par0259+0032_194, and
show the photometric redshift fits when the $J_{125}$ photometry is excluded
(solid green curve) and when a lower limit is asssumed as described in the
text (dashed black curve).  Figure 11: The candidate plots for the remaining 7
candidates. Plots are displayed as in Figure 10.

In Figures 10 and 11 we show the best-fitting SED and photometric redshift
PDFs for the HST photometry (purple dashed) and the HST+IRAC photometry when
available (blue solid). We also show the best-fitting template obtained by
limiting EAZY to $z<7$ as a red dash-dotted curve. The lower-redshift ($z<7$)
templates are ruled out by the upper limits measured blueward of the expected
Lyman break, especially those in the $Y$-bands for candidates at $z<10.5$ and
in the $J$-bands for candidates at $z>10.5$. The addition of the Spitzer/IRAC
photometry also significantly helps, especially the upper limits measured in
the 3.6µm images for candidates at $z>10$. The IRAC non-detections rule out
dusty galaxies with redder $H$-[3.6µm] colors in favor of higher-redshift
sources with a flat spectrum or bluer color. Postage stamps of all available
filters are shown along the top, where we display the Galfit residual maps for
each individual IRAC observation. All modeled fluxes are removed from the
residual maps except that attributed to the candidate.

In the following subsections, we briefly discuss each candidate.

####  4.1.1 Par1033+5051_116

Candidate Par1033+5051_116 ($m_{160}=25.99$) is detected in both $J_{125}$ and
$H_{160}$. As a two-band detection, we consider this candidate to be higher-
confidence than the single-band detections in the sample, both because it is
less likely to be a spurious source and because its photometric redshift is
better constrained. With a photometric redshift of $z_{\mathrm{phot}}=8.28$,
this candidate is at the low end of our sample, yet the 68% interval extends
out to $z\sim 8.7$.

####  4.1.2 Par0440-5244_497 & 593

Candidate Par0440-5244_497 ($m_{160}=25.47$) is also a two-band detection,
with a $z_{\mathrm{phot}}=8.53$. A second candidate in our sample,
Par0440-5244_593 ($m_{160}=26.46$) is detected in the same HST pointing at a
distance of $\sim$1$\farcm$25 from source 497, corresponding to $\sim$0.35
pMpc at $z=8.5$. However, as Par0440-5244_593 is a single-band detection with
a photometric redshift of $z_{\mathrm{phot}}=11.07$, the proper distance
between the two sources is $\sim$58 pMpc. Here we have adopted each of their
photometric redshifts and used the average $z_{\mathrm{phot}}$ to convert
their comoving distance to a proper distance. Given the redshift separation,
we do not expect these sources to be associated with the same overdensity (see
discussion in Section 7.2).

####  4.1.3 Par0713+7405_95

The third candidate in Table 5, Par0713+7405_95 ($m_{160}=25.74$), is also a
two band detection. While the photometry in Figure 10 shows a detection in
$I_{600}$, we note that it is not significant with S/N $<1.5$. The inclusion
of the IRAC photometry helps to narrow the width of the primary redshift PDF
peak, placing $z_{\mathrm{phot}}=9.23$. However, the rejection of the lower-
redshift solution is driven by the non-detection in $Y_{105}$.

####  4.1.4 Par0456-2203_473

The IRAC deblending for Par0456-2203_473 is uncertain due to its close, bright
neighbor. The flux remaining in the IRAC residual maps has been attributed to
the candidate, but it is not fully centered in the blue box that indicates the
size of the HST stamp. This flux may therefore be partially or fully coming
from the neighboring source. With such a close neighbor, any measured fluxes
at IRAC wavelengths will be uncertain until the much higher resolution
JWST/NIRCam can be used to separate the light from these sources. However, the
HST-only photometry yields a high-redshift solution that is still preferred
even when these relatively bright IRAC fluxes are included. We note that
Par0456-2203_473 has a bright $H_{160}$ magnitude ($m_{H160}=24.52$), which
could make it one of the brightest candidates detected at $z>8$. We calculate
an absolute UV magnitude for this source of $M_{1500}=-22.86$, $\sim$0.2
brighter than the $z\sim 10$ candidate presented by Morishita et al. (2018) in
the BoRG[$z$910] fields. However, as discussed in Section 6.1, this source is
strongly lensed by its very close neighbor, and is likely much fainter
intrinsically. This source is also the only one in our sample that received
the aperture correction for stretched apertures described in Section 2.4.2.
The $H_{160}$ magnitude in Table 5 is an aperture-corrected flux calculated in
a circular aperture, rather then the elliptical Kron apertures used for the
other candidates. We note that we measure a large half-light radius for this
candidate, $r_{1/2}=0\farcs 33$. This $r_{1/2}$ is consistent with the radius
measured by Holwerda et al. (2015) for low-redshift interlopers, yet is not
inconsistent with a high redshift nature (Holwerda et al., 2020). Also, as a
strongly lensed source, this 0$\farcs$33 radius may not represent the
intrinsic size of this candidate.

####  4.1.5 Par0259+0032_194

A fifth candidate, Par0259+0032_194 ($m_{160}=24.79$), is potentially another
two-band detection, though there is a cluster of bad pixels in the $J_{125}$
image at the position of the source (there is no such cluster affecting
$H_{160}$). The $J_{125}$ flux of this source was masked out in our
photometric catalog, and the source is treated as having no $J_{125}$ coverage
in our analysis (sample selection, photometric redshift fitting, comparison
with SpeX spectra). However, the removal of $J_{125}$ leads to a poorly-
constrained redshift solution with the source placed at a higher redshift than
it likely should, considering it may very well be detected in $J_{125}$. We
re-measure the $J_{125}$ flux of Par0259+0032_194 after masking out the bad
pixels (using Source Extractor in dual image mode with $H_{160}$ as before).
We use this $J_{125}$ flux to measure a new photometric redshift,
$z_{\mathrm{phot}}=9.36$, which is quoted in Table 5. We adopt this redshift
PDF (grey “$H_{160}$ and $J_{125}$” curve in Figure 10) in estimating the
candidate’s magnification and absolute magnitude posterior (Section 6.1) and
constructing the luminosity function (Section 6.4). We note that this
$J_{125}$ flux is likely a lower limit, and a brighter $J_{125}$ flux would
indicate a stronger spectral break and therefore more strongly reject the
smoother slope of a lower-redshift, red galaxy. In Figure 10, we show the
original SED and redshift PDF for Par0259+0032_194 (i.e., excluding the
$J_{125}$ flux) in green and those including the estimated $J_{125}$ lower
limit in black, using a different color scheme as the rest of the figure to
signify this special case.

####  4.1.6 Par335_251

Included in our sample is the first $z>8$ candidate from the WISP survey,
Par335_251 ($m_{160}=24.74$). It has a photometric redshift of
$z_{\mathrm{phot}}=10.53$ and is detected only in the $H_{160}$, as expected
for WISP fields (see discussion in Section 6.3). As can be seen in Figure 10,
it lies on the UVIS chip gap in $V_{606}$, and so we treat the source as
having no $V_{606}$ coverage. As Par335 has IRAC 3.6µm coverage, we have this
additional constraint for the WISP source, and so we include it in our sample.
Additionally, as a slitless spectroscopic survey, WISP observed Par335 with
both of WFC3’s IR grisms: G102 ($8.8-1.1$ µm, $R\sim 201$) and G141
($1.07-1.7$ µm, $R\sim 130$). The integration times were tuned to achieve
approximately uniform sensitivity at all wavelengths for a line of a given
flux, and for this field are $\sim$8800 seconds in G102 and $\sim$3800 seconds
in G141.

Figure 12: The one-dimensional WISP spectra in the WFC3 grisms G102 (blue) and
G141 (orange) of candidate Par335_251. The light shaded regions show the
1$\sigma$ uncertainties. The best-fitting high-redshift solution
($z_{\mathrm{phot}}=10.53$) would place Ly$\alpha$ at
$\lambda_{\mathrm{obs}}\sim 1.4$µm (blue dashed vertical line). The best-
fitting low-redshift solution would place [O II] at
$\lambda_{\mathrm{obs}}\sim 1.32$µm (red dash-dotted line). We show a
horizontal rectangle for each emission line weighted by the redshift PDFs and
extending to their 95% intervals. There are potential emission lines at $\sim$
12200Å (S/N=2.3) and 14000Å (S/N=1.5) that are broadly consistent with the
$z$$\sim$2 redshift PDF, though this low-redshift solution is disfavored with
a $\Delta\chi^{2}=20.41$. The fits to each line are described in the text and
displayed in the top panels. The 14000Å feature would be consistent with
Ly$\alpha$ at $z=10.53$, though a Ly$\alpha$ detection at this redshift is
unlikely due to the relatively short integration times. While there is no
strong indication that this candidate is at a lower redshift, its true nature
is not clear from this spectrum.

While these spectra will be far too shallow to detect an emission line or
continuum for a source at $z\sim 10.5$, it is instructive to explore them for
any spectral features that may indicate the source is at a lower redshift. We
downloaded the extracted, calibrated spectra from the WISP page on MAST (see
footnote 4) and plot them in Figure 12. Assuming the best-fitting redshifts
for the high- and low-redshift solutions, we show the expected location of
Ly$\alpha$ at $z=10.53$ and [O II]$\lambda$3727 at $z=2.54$. There are two
possible emission lines at $\lambda_{\mathrm{obs}}\sim 12200$Å and
$\lambda_{\mathrm{obs}}\sim 14000$Å that lie on either edge of the expected
range for [O II] given the low redshift solution.

We measure these lines using a Monte Carlo analysis, assigning fluxes at each
wavelength step that are pulled randomly from a normal distribution centered
at the flux value and with $\sigma$ equal to the flux uncertainty. We then fit
the fluxes with a Gaussian profile using least-squares minimization to
determine the best fit, and repeat this process 1000 times. For each line, the
median of all 1000 fits is displayed in the top panels of Figure 12. We define
the S/N of the line as the flux of the median fit divided by the standard
deviation of the fits. The potential emission lines at 12200Å and 14000Å have
S/N=2.3 and 1.5, respectively. We note that the sky background has been over-
subtracted in the spectrum, and so we fit a new continuum (black horizontal
dashed line) from which we measure the strength of each line. If fit instead
with a continuum at flux=0, the S/N of each line would decrease.

Both lines have S/N$<$3, and lie far from the expected observed wavelength for
[O II]$\lambda$3727 given the low redshift solution’s narrow PDF.
Additionally, the $\Delta\chi^{2}$ between the best-fit solution at $z\sim
10.5$ and the minimum $\chi^{2}$ at lower redshift is 20.41, indicating that
the low-redshift solution is disfavored. The potential line at
$\lambda_{\mathrm{obs}}\sim 14000$Å would be very consistent with Ly$\alpha$
at a redshift of $z=10.53$. However, such a detection would require hours more
integration time, as in Zitrin et al. (2015) and Larson et al. (2022). While
the WISP spectrum does not exhibit clear signs of being at a lower redshift,
it cannot help constrain the nature of Par335_251.

####  4.1.7 Remaining Single-Band Detections

The remaining five candidates:

  * •

Par0956-0450_684 ($m_{160}=24.88$, $z_{\mathrm{phot}}=10.54$)

  * •

Par0843+4114_120 ($m_{160}=26.39$, $z_{\mathrm{phot}}=10.98$)

  * •

Par1301+0000_37 ($m_{160}=25.42$, $z_{\mathrm{phot}}=11.02$)

  * •

Par0756+3043_92 ($m_{160}=26.03$, $z_{\mathrm{phot}}=11.03$)

  * •

Par0926+4536_155 ($m_{160}=26.37$, $z_{\mathrm{phot}}=11.15$)

are all single-band detections in the HST imaging. All five have IRAC coverage
at both 3.6µm and 4.6µm, and are non-detections in both bands. While the non-
detections in the bluer HST filters as well as both IRAC channels help to rule
out lower-redshift contaminants, there is an increased concern that these
single-band detections are spurious sources (e.g., hot pixels or cases of
persistence). Three of these candidates are from the HIPPIES GO 12286 dataset,
which has been noted to suffer from increased levels of persistence and
contamination (e.g., Morishita, 2021). We have done our best to identify and
remove spurious sources (see Sections 3.4 and 3.3), yet the possibility
remains that some or all of these candidates are contaminants in our sample.
Spectroscopic confirmation or imaging in more filters with NIRCam is needed to
confidently evaluate the high-redshift nature of these sources.

###  4.2 Comparison with Previous Studies in the Same Fields

The datasets we include in our analysis are from surveys that were designed to
detect high-redshift galaxies. The BoRG survey in particular was optimized for
searches for $z\gtrsim 8$ sources, and several teams have previously reported
results for galaxies in these fields. In this section, we discuss results from
the previous studies that were performed in the same HST fields that we
consider in this paper. As this is the first search for $z>8$ candidates in
HIPPIES GO 12286 or WISP fields, we focus here on results in the BoRG[$z$8]
fields.

We begin by discussing the BoRG[$z$8] results aimed at searching for galaxies
at $z\sim 8$ as $Y_{098}$ dropouts (Trenti et al., 2011; Bradley et al., 2012;
Schmidt et al., 2014a). These studies require a significant detection in
$J_{125}$ and use $H_{160}$ as a second detection band to rule out spurious
sources. We are focusing our search at higher redshifts and allow for single-
band detections, and so do not expect to recover many of the candidates
presented in these papers. However, we do find two candidates in common, which
is unsurprising because the redshift PDFs of some sources in our sample extend
to $z\sim 8$. Our two lowest-redshift candidates, Par1033+5051_116 and
Par0440-5244_497, were previously published by Bradley et al. (2012). For both
candidates, we measure slightly brighter $H_{160}$ magnitudes:
$m_{H160}=25.99$ (Par1033+5051_116) and 25.47 (Par0440-5244_497) compared with
26.2 and 25.9, respectively, based on the $J_{125}$-$H_{160}$ colors reported
in Bradley et al. (2012). While Bradley et al. (2012) do not calculate
photometric redshifts for their candidates, the completeness they present for
a representative field peaks at $z\sim 7.8-8$, which is consistent with the
68% intervals we measure for our redshift PDFs. We do not recover the
remainder of the candidates presented in these three papers because our
selection function is more sensitive to galaxies at higher redshifts than
those identified as $Y_{098}$-dropouts. We also focus on a slightly brighter
magnitude range than these BoRG[$z$8] studies, requiring $m_{H160}\leq 26.5$
while many of the candidates presented in these papers are in the range
$26.7\lesssim m_{H160}\lesssim 27.2$.

Next, Morishita et al. (2020) perform a search for quasars at $z\sim 8$ in the
BoRG[$z$8] fields (among others), focusing on the identification of Lyman
break candidates that have point-source morphologies rather than resolved
shapes. Morishita et al. (2020) also employ a combination of color cuts and
detection significance criteria in their sample selection, but additionally
require that the $\geq$70% of the source redshift PDFs be at $z>6.5$.
Morishita et al. (2020) are again focusing on a slightly lower redshift range
than we are, and so it is not surprising that our samples do not overlap. We
do not select their one candidate in a BoRG[$z8$] field because it is detected
in $Y_{105}$ with a S/N$>2$. We do find one candidate in the same field as
they do (Par0926+4536_155), but our candidate is $\sim$2′ away from the one
they identify.

While the BoRG[$z$8] survey was designed to detect galaxies at $z\sim 8$ as
$Y_{098}$-dropouts, Bernard et al. (2016, hereafter B16) performed a search
comparable to ours for $z\sim 9-10$ galaxies in these fields. They focused on
the 62 fields with an exposure time in $H_{160}$ of at least 900 seconds, a
survey area covering 293 arcmin2. (We include 66 pointings and make no
exposure time cut.) B16 identified three $z\sim 10$ candidates. They perform
their search using the Lyman break technique to identify $J_{125}$ dropouts,
implementing a $J_{125}-H_{160}>1.5$ color selection that will select a strong
$J_{125}$$-$$H_{160}$ break while still allowing for $J_{125}$ detections.
This color selection is coupled with a detection requirement of
S/N${}_{H160}\geq 8$ and non-detections required in the $V$ and $Y_{098}$
filters (S/NV and S/N${}_{Y098}<1.5$). They also implemented a threshold of
CLASS_STAR$<0.95$ and performed a visual inspection of their candidates. This
source selection method is in contrast to our process, which we base on
photometric redshifts rather than color-cuts (Section 3). As discussed in
Section 2.5, we also find that the accuracy of the CLASS_STAR flag decreases
with magnitude, and so we did not incorporate it into our selection but
instead compare the colors of our sources to identify stellar contaminants. We
also note that B16 did not have access to Spitzer/IRAC imaging of their
candidates, but applied for IRAC follow-up which we are able to include in our
analysis here (Spitzer PID 12058, PI Bouwens).

B16 identified a further three candidates as contaminants based on source
size. Following Holwerda et al. (2015), the authors adopt 0$\farcs$3 as an
upper limit for the half-light radius of their candidates and reject larger
sources as low-redshift interlopers. On the other hand, Holwerda et al. (2020)
find that while such a size criterion could remove up to 75% of low redshift
interlopers, it may also reject $\sim$20% of real, bright galaxies in the
targeted redshift range. Source size and morphology is another area where the
higher resolution of the JWST/NIRCam imaging will significantly improve our
understanding of the size-luminosity relation out to $z\gtrsim 10$. We do not
apply an upper limit size cut to our sample selection, though our candidates
are fairly compact with (non-PSF-corrected) half-light radii in the range
0$\farcs$11$-$0$\farcs$19 (with the exception of Par0456-2203_473, see Section
4.1.4).

We do not recover any of the three candidates presented in B16. Two of their
candidates (borg_0240-1857_25 and borg_0456-2203_1091) satisfied our selection
criteria but were rejected during visual inspection in category 1e (a hot
pixel, cosmic ray, or source with a strange morphology). We identify these two
sources in Figure 19. Our rejection of these two candidates illustrates the
necessarily subjective nature of high-redshift candidate searches with HST,
where potentially real candidates may be removed in an effort to clean samples
of spurious sources. Interestingly, we identify a different candidate in one
of the same fields, Par0456-2203_473. Our candidate is located $\sim$56″away
from borg_0456-2203_1091, yet was not selected by B16. We did not detect the
third candidate presented in B16 (borg_1153+0056_514) with enough significance
in $H_{160}$ to meet our S/N${}_{H160}>5$ criterion (S/N=4.78) and so this
source was not selected in our sample. The Source Extractor parameters used to
detect objects in an image, the choice of aperture size and shape, and the way
in which the noise is calculated can all affect the detection significance.
These discrepencies between the sample we present in Section 4.1 and that of
B16 highlight how dependent $z>8$ candidate samples are on the sample
selection methods used.

##  5 Low-Redshift Contamination

The datasets we consider in this paper are relatively shallow compared to
those of the well-studied fields such as CANDELS or the Hubble Frontier Fields
that are the best current options for detecting $\leq$$L*$ high-redshift
galaxies. While these parallel datasets are particularly well-suited to
detecting rare bright galaxies needed to explore the bright end of the
luminosity function, they are likely not deep enough to detect the far more
common faint galaxies at lower-redshifts. The possibility remains that our
selection criteria are in fact identifying faint, low redshift galaxies.
Faint, red galaxies at $z\sim 2-3$ are of particular concern. Whether passive
or star-forming and dusty, their steepening red SEDs may only rise above our
detection thresholds in the reddest filters ($J$ and $H$). Extreme emission
line galaxies at lower redshift are another possible source of contamination,
where high equivalent width lines such as [O II], [O III], or H$\alpha$ could
boost the $H_{160}$ flux at $z\sim 3$, $z\sim 2$, or $z\sim 1.3$, respectively
(e.g., Atek et al., 2011; Faisst et al., 2016). In both cases, photometry in
IRAC bands can help identify these sources as their SEDs likely increase past
2µm. However, a faint contaminant with an emission-line-boosted $H_{160}$ flux
may fall below the IRAC detection limit. Additionally, one of our candidates
(Par0259+0032_194) does not have IRAC coverage, and we consider the IRAC
deblending for another candidate (Par0456-2203_473) to be unreliable. While
Spitzer/IRAC imaging should help break these redshift degeneracies, the
shallow depths of the IRAC imaging as well as the uncertainties in the IRAC
deblending complicate the picture. Our conservative selection criteria and the
imaging currently available at $\sim 3-5$µm are not sufficient to conclusively
rule out such contaminants. In the following sections, we therefore explore
the low-redshift contamination fraction in our sample from faint galaxies
falling below our detection thresholds.

###  5.1 Estimate from Stacking

Figure 13: Stacks of image postage stamps and the photometric redshift
probability distributions of all candidates in our sample. We separate the
stamps by dataset in order to combine images with the same pixelscales,
showing BoRG[$z$8] in the top row and HIPPIES GO 12286 in the second row.
These stamps are 5″ on a side, displayed on a zscale interval with no image
smoothing. The black circle of $r=0\farcs 5$ indicates the expected position
of source flux. In the bottom row we show the stacked redshift PDF of all
sources in our sample (left) and a stack of all $V$, $I$, and $Y$ imaging for
all candidates (including WISP) resampled onto a common pixel scale of
0$\farcs$01/pixel (right). In all cases, including the resampled stack, no
source flux is visible in the optical and $Y$ filters at the expected
position, indicating that our sample is not dominated by lower-redshift
contaminants.

The optical imaging in these HST datasets may be too shallow to detect flux
from individual low-redshift sources, and so we stack postage stamps of the
candidates in order to check for any optical detections. We create several
stacks, separating the candidates by dataset to avoid introducing any
artificial features by combining images on different pixel scales. We also
separate the imaging by filter, creating one stack for each dataset that
includes all available images obtained with the WFC3/UVIS and ACS/WFC cameras
($V_{606}$ and $I_{600}$ for BoRG[$z$8], and $I_{600}$ for HIPPIES GO 12286),
a stack of all WFC3/IR $Y$ images as an additional dropout band, and stacks of
all $J$ and $H$ images as the candidate detection bands. We do not include the
candidate from the WISP survey in these stacks because it is the only one on
the 0$\farcs$13 pixel scale. In each case, the images are combined using a
weighted mean in order to down-weight any bad pixels. The image stacks are
displayed in the top two rows of Figure 13, where we indicate the number of
images in each filter that contributed to the stack. As can be seen in the
left two panels of each row, there is no appreciable source flux at the
expected position in the dropout bands.

We next estimate the detection significance at the center of each stack using
$r=0\farcs 2$ circular apertures. We place a grid of unique apertures across
each stamp, using a segmentation map created by combining the individual maps
for each input field to avoid source flux. The $1\sigma$ noise in the stamp is
then taken as the standard deviation of a Gaussian distribution fit to the
distribution of background aperture fluxes. With the aperture fluxs measured
at the center of each stamp, we find a signal to noise of $\mathrm{S/N}<1$ in
both sets of $V+I$ and $Y$ stamps, confirming that there is no significant
optical detection in the stacked imaging. These imaging stacks provide at best
an upper limit estimate of the low redshift contamination in the sample,
suggesting that the majority of our sample is not comprised of contaminants.
For reference, the S/N in $H$ ($J$) is 24.1 (3.7) and 65.5 (9.3) for the
BoRG[$z$8] and HIPPIES GO 12286 stacks, respectively.

As an extra check, we resample all available $V$, $I$, and $Y$ imaging,
including the $I_{814}$ imaging from the WISP candidate, to a much finer pixel
scale of 0$\farcs$01/pixel. This resampling then allows us to combine all 25
images into a single, summed stack, which is shown in the bottom right corner
of Figure 13. If a majority of the candidates in our sample were in fact
lower-redshift interlopers, we would expect to detect flux in this combined
stack of 25 images. However, there is no visible signal in the black circle in
the bottom right panel of Figure 13, indicating that lower-redshift
contamination does not dominate our sample.

In the bottom row of Figure 13, we also show the stacked redshift PDFs for all
13 sources in our sample. There is a non-zero portion of the PDF at $z\sim
2-4$, with a peak around $z\sim 2.5$. This low-redshift peak is consistent
with the detection of the 4000Å break in lieu of the Ly$\alpha$ break between
$J_{125}$ and $H_{160}$. Alternatively, if high equivalent width emission line
galaxies are contributing to that low-redshift peak, it could correspond to
strong [O III] in $H_{160}$ up to $z\sim 2.4$, with [O II] transitioning into
the $H_{160}$ filter at $z\sim 2.75$.

The portion of the stacked PDF that lies at $z<7$ amounts to only 1.6% of the
total integrated PDF. However, as we have not used luminosity priors in our
EAZY runs (see Section 3.1), the redshift PDFs do not take into account the
relative abundances of high- and low-redshift galaxies in the universe.
Therefore, we cannot use the redshift PDF to estimate the contamination from
lower-redshift galaxies in our sample. Instead, in the next section we aim to
quantify the contamination fraction by artifically dimming low-redshift
sources and attempting to reselect them as high-redshift candidates.

###  5.2 Catalog-Level Estimate of Faint Lower-redshift Interlopers

We begin by identifying a high-fidelity sample of galaxies at low redshift.
From the full Source Extractor catalog, we select all sources with
S/N${}_{H}>5$ and $H_{160}<23.5$ mag. After running EAZY on this sample of
bright sources, we select all those with 70% of their redshift probability
distribution contained in the range $1\leq z\leq 4$. We consider only resolved
sources in an attempt to avoid including in this analysis any stars that have
been incorrectly assigned $z>0$ photometric redshifts. Using the size
measurements presented in Section 3.7, we select sources with a half-light
radius larger than the median $r_{1/2}+1\sigma$ measured for stars in the
corresponding dataset. There are a total of 622, 142, and 335 low-redshift
sources – those that are bright, resolved, and reliably in the range $1\leq
z\leq 4$ – in the BoRG[$z$8], HIPPIES GO 12286, and WISP datasets,