respectively.

For each HST field containing a candidate high-redshift galaxy in our sample,
we create a catalog of dimmed, low-redshift galaxies by assigning a new
$H_{160}$ magnitude to each real source. The dimmed $H_{160}$ magnitudes are
pulled randomly from a truncated Gaussian distribution with $\mu=26.5$ and a
maximum at $m=27$, 0.5 magnitudes fainter than the 5$\sigma$ requirement in
our selection criteria. We dim the fluxes in all filters by the same value,
$\Delta m=H_{\mathrm{obs}}-H_{\mathrm{dimmed}}$. This process maintains the
colors of a $z\sim 2$ galaxy while moving all photometry to fainter
magnitudes. We create multiple realizations of each input low-redshift source,
resulting in a catalog of 5000 dimmed sources for each field. We measure the
median flux uncertainties in magnitude bins for all sources in the Source
Extractor catalog that are detected in the given HST field, and perturb the
dimmed fluxes by the error for the appropriate magnitude bin. In this way, we
incorporate field-to-field variations in depths while accounting for
photometric scatter as a function of magnitude. Finally, we run EAZY on the
new photometry and reselect sources as those with $H_{160}\leq 26.5$ and 70%
of their redshift PDF at $z>8$.

Figure 14: An example of a low-redshift source from Par0259+0032 that is
selected as a high-redshift candidate after the fluxes are dimmed as described
in the text. In the large panel on the left, we plot the original and dimmed
photometry as orange circles and red squares, respectively. The best-fitting
EAZY templates are also displayed, with the corresponding $\chi^{2}$
distributions in the top right panel and the redshift PDFs in the bottom right
panel. We show the original (undimmed) HST imaging in 5″stamps along the top
left. Dimming the fluxes in all bands preserves the source’s colors, but
results in the source dropping out of the bluer bands, and the new photometry
is better fit at $z\sim 9$. Through this analysis, we expect low contamination
rates ($<$1%-12%) from faint lower-redshift galaxies in our sample.

For each field, we then look at the total number of the 5000 sources that were
“recovered” as high-redshift candidates. There were no dimmed galaxies
selected as high-redshift sources in five of the 12 fields that contain our
candidates (Par0440-5244, Par0713+7405, Par0756+3043, Par1033+5051 and
Par1301+0000), implying a near-zero contamination rate. For the remaining
fields, we select between $1-4$ dimmed sources. Figure 14 shows an example of
one such low-redshift source, with the original (undimmed) HST image stamps of
the source as well as the original and dimmed photometry plotted as orange
squares and red circles, respectively. We also show both sets of EAZY template
fits, redshift PDFs, and $\chi^{2}$ distributions. As can be seen here, the
dimming process maintains the $J-H$ color of the source within the photometric
scatter, but the bluer fluxes have dropped out of the optical bands.

We use equation 4 from Finkelstein et al. (2015b) to estimate the
contamination fractions in each field that correspond to the number of
recovered sources:

$F=\frac{N_{\mathrm{dim,sel}}}{N_{\mathrm{dim,low}z}}\frac{N_{\mathrm{tot,low}z}}{N_{\mathrm{high}z}}.$
(2)

Here, $N_{\mathrm{dim,sel}}$ is the number of dimmed low-redshift sources that
were selected as high-redshift galaxies, and $N_{\mathrm{dim,low}z}$ is the
number of input dimmed sources (5000 in all fields). The ratio
$N_{\mathrm{tot,low}z}/N_{\mathrm{high}z}$ accounts for the relative number
densities of low- and high-redshift sources in our survey volume, where
$N_{\mathrm{high}z}$ is the number of high-redshift candidates detected in the
dataset (i.e., 7 candidates in BoRG[$z$8] fields, 5 in HIPPIES GO 12286, and 1
in WISP), and $N_{\mathrm{tot,low}z}$ is the total number of faint, low-
redshift galaxies identified in the catalog of each dataset. For the latter
number, we select low-redshift sources with 70% of their integrated redshift
PDF contained in the range $1\leq z\leq 4$ and magnitudes in the range
$24<m_{H160}\leq 27$. We include sources below the threshold of $m_{H160}\geq
26.5$ in order to account for the prevalence of fainter low-redshift galaxies
that may be photometrically scattered into our sample magnitude range. In an
effort to avoid stars and spurious detections and to ensure quality
photometric redshift fits, we also require a half-light radius larger than the
median measured for stars in the catalog as well as a $\geq$5$\sigma$
detection in $H_{160}$. We note that these selections, and especially the last
two, are likely to miss real faint, low-redshift sources in our catalog.
However, we find that without these requirements, we incorrectly include many
stars and marginally-detected sources with poorly-constrained redshift PDFs in
our total for $N_{\mathrm{tot,low}z}$.

As an example, we find that BoRG[$z$8] field Par0456-2203 has:
$N_{\mathrm{dim,sel}}=1$ with $N_{\mathrm{tot,low}z}=3342$ in the BoRG[$z$8]
catalog, and $N_{\mathrm{high}z}=7$, resulting in a contamination fraction of
$F=0.095$. We find contamination rates of 1.4% (Par0926+4536), 2.2%
(Par0956-0450 and Par0259+0032), 2.9% (Par0843+4114), 9.5% (Par0926+4000 and
Par0456-2203), and 12.7% (Par335). The fields with $N_{\mathrm{dim,sel}}=0$
have contamination rates of $<$1% (Par0713+7405) and $<$9.5% (Par1033+5051,
Par0440-5244, Par0756+3043, and Par1301+0000). We again note that these
contamination rates may be underestimated. We have limited our analysis to
galaxies in the range $1<z<4$ because the majority of the low-redshift
portions of our candidates’ photometric redshift PDFs is in this range. Yet we
are therefore not accounting for low-redshift contamination from sources
outside of this range. Also, the SEDs of sources in the catalog are
constrained by only a handful of HST photometric measurements, resulting in
uncertainties in the photometric redshifts measured for the bright, low-
redshift sources that are chosen to be dimmed, the dimmed sources that are
then ‘recovered’ as high-redshift galaxies ($N_{\mathrm{dim,sel}}$), and the
faint, low-redshift sources ($N_{\mathrm{tot,low}z}$). However, with at most
four bands of photometry, this analysis represents our best attempt at
identifying contamination from faint, low-redshift galaxies.

As the estimated contamination rates are small, and as we expect the
contamination rates measured in the bins of magnitude we use when calculating
the UV luminosity function to be even lower, we do not reduce our observed
number densities by the fractions reported here.

##  6 The Bright End of the Luminosity Function at $z\sim 9$

In this section we calculate the rest-UV luminosity function at $8.5\lesssim
z\lesssim 11$. This measurement first requires correcting the candidate
luminosities for any magnification imparted by nearby neighbors (Section 6.1),
and estimating the completeness of our sample as a function of magnitude
(Sections 6.2 and 6.3).

###  6.1 Source Magnification

Gravitational lensing can have the effect of boosting the brightness of
sources such that they are detected in a flux-limited sample when
intrinsically they would be too faint. Lensing can benefit high-redshift
surveys by bringing the fainter portions of the galaxy population within reach
of HST surveys (e.g., Ellis et al., 2001; Coe et al., 2013, 2015; Schmidt et
al., 2014b; Atek et al., 2015; Bouwens et al., 2017; Livermore et al., 2017),
providing more constraining power for the faint end slope of the luminosity
function to higher redshifts than would otherwise be possible. However, the
magnification of sources can also significantly change the shape of the
measured luminosity function, especially at the bright end where the low
number counts are easily boosted by an increase in observed luminosity. Many
of the candidates in our sample have close neighbors, and five in particular
have one or more bright neighbors within $\sim$3″. We now estimate the
magnification that each candidate experiences, following the steps laid out by
Mason et al. (2015b) for strong lensing.

The observed magnitude of a candidate will be increased by a factor
$-2.5log_{10}(\mu)$ by a bright, close neighbor. The magnification, $\mu$,
depends on the redshifts of the candidate and the lensing source, the velocity
dispersion of the lens, and the projected distance between the lens and the
candidate in the image plane. Following Mason et al. (2015b), we assume each
of the potential lensing sources can be modeled as a singular isothermal
sphere, where the Einstein radius of the lens is proportional to the square of
the velocity dispersion of the lensing galaxy ($\sigma^{2}$) (see equation 4
of Mason et al., 2015b).

We estimate $\sigma$ using the empirical, redshift-dependent Faber-Jackson
relationship that Mason et al. (2015b) dervive using BoRG[$z$8] fields (see
their Table 1). For this estimate, we use the EAZY-calculated photometric
redshifts for the neighboring sources, including the IRAC photometry where
available to improve on the photometric redshift constraints. We obtain IRAC
photometry for each neighbor as described in Section 3.5, adopting the Galfit
model magnitudes for detections and using $r=2.5$ pixel aperture fluxes for
non-detections. When calculating $\theta_{\mathrm{ER}}$, ideally we would use
the apparent magnitudes from the same filters as Mason et al. (2015b) for each
redshift bin. However, we find that the majority of the neighboring sources
have a significant detection (S/N$\geq$3) only in $H_{160}$, and so we use
$H_{160}$ apparent magnitudes for all neighbors, regardless of redshift. As
can be seen in Figures 10 and 11, some of the candidates have multiple bright,
close neighbors. We assume that each neighboring source is an independent
lens, such that the total magnification $\mu$ for a candidate is the product
of the magnification from each neighboring source, $\mu_{i}$. This method
assumes that the neighbors each lens the light from the candidate along the
line of sight independently, an assumption that we discuss in more detail
below.

As in Rojas-Ruiz et al. (2020), we use a Monte Carlo approach to calculate
$\mu$, the corrected $H_{160}$ apparent magnitude, and the absolute rest-UV
magnitude, $M_{\mathrm{UV}}$, of each candidate. We perform 5000 calculations
of $\theta_{\mathrm{ER}}$ and $\mu_{i}$ for each neighbor within 5″ of a
candidate, increasing the number of realizations from 1000 in Rojas-Ruiz et
al. (2020) because of the broader redshift PDFs in our catalog. For each
realization, we pull the redshift for the candidate and each neighbor randomly
from their respective redshift PDFs. We randomly draw the $H_{160}$ fluxes for
the candidate and neighbors from Gaussian distributions centered on the
observed value and with a standard deviation equal to the flux uncertainties.
In this way, we are including the uncertainties in the photometric redshifts
and the flux uncertainties in our calculations of $\mu$ and $M_{\mathrm{UV}}$.
For the rest-UV magnitude, we use the fact that a rest-frame 1500Å wavelength
falls near the middle of the $H_{160}$ filter in the $z\sim 9.5$ observed
frame.

Finkelstein et al. (2022) show that including all nearby neighbors in these
calculation biases the resulting magnification towards higher values of $\mu$.
Specifically, they show that on average this approach results in $\mu\sim 1.4$
along randomly-selected lines of sight where the expected magnification is
$\mu=1$. Following their example, we therefore only consider the magnifaction
contribution from neighbors with a median $\mu_{i}>1.4$. The $\mu_{i}$ of any
neighbor with $\mu_{i,\mathrm{median}}<1.4$ is set to unity. As a result of
our Monte Carlo process, we have a distribution of 5000 $\mu$’s, corrected
$H_{160}$ fluxes, and $M_{\mathrm{UV}}$ values. We take the median of each
distribution as the parameter estimate, and adopt the 68% interval as the
uncertainty.

We find that three candidates experience a significant magnification:
Par0713+7405_95 with $\mu=1.45\pm 0.13$, Par0259+0032_194 with $\mu=3.81\pm
11.8$, and Par0456-2203_473 with $\mu=15.68\pm 12.50$ (note that this
candidate also has the largest observed $r_{1/2}$). This intermediate-to-
strong lensing is unsurpising as each of these three candidates has a bright,
close neighbor. However, we note that the magnifications of the two candidates
with the strongest lensing are highly uncertain, with $\sigma_{\mu}>\mu$ for
Par0259+0032_194. We report the median $M_{\mathrm{UV}}$ for each candidate in
Table 5, adopting the magnification-corrected magnitude for the three lensed
candidates and providing the 68% intervals that account for the uncertainties
in photometric redshift, $H_{160}$ fluxes, and magnification where applicable.
We use the distribution of $M_{\mathrm{UV}}$ values in calculating the
luminosity function in Section 6.3.

###  6.2 Simulated Sources

In order to quantify the incompleteness in our sample of high-redshift
galaxies caused by our selection criteria, we create a catalog of simulated
sources and add them to the real images from each WFC3 pointing. We then
perform the same steps of source detection, photometry, and sample selection
described in Sections 2.4 and 3.2, treating the images containing simulated
sources identically to how we treat the original images.

We create a catalog of $10^{5}$ sources, with each source assigned a redshift,
spectral template, observed flux, size, and shape from parameter distributions
as follows. The source redshifts are pulled randomly from a uniform
distribution in the range $7\leq z\leq 12$. Each source is then assigned a
spectral template from the Bruzual & Charlot (2003) library with either a
Salpeter (1955) or Chabrier (2003) initial mass function, and an
exponentially-declining star formation history with a characteristic timescale
of either $\tau=0.01$, 0.5, or 5.0 Gyr. These timescales are adopted to
approximate a single burst of star formation at one extreme ($\tau=0.01$ Gyr)
and a constant star formation history at the other ($\tau=5$ Gyr). The
template metallicities are one of 0.02, 0.2, 0.4, or $1.0Z_{\odot}$, and are
pulled randomly from a lognormal distribution that peaks at $0.2Z_{\odot}$.
The ages of the stellar populations are similarly taken from a tweaked
lognormal distribution that peaks at 10 Myr and decreases with an enchanced
tail towards older ages to ensure that some older ($>400$ Myr) populations
exist. The ages available for each simulated source are restricted to be
between 10 Myr and the age of the universe at the given redshift. To each
spectral template, we then apply dust extinction due to the interstellar
medium using the Calzetti et al. (2000) extinction law and a color excess
$E(B-V)$ pulled randomly from a truncated Gaussian distribution with
$\mu=0.1$, $\sigma=0.15$, $E(B-V)_{\mathrm{min}}=0$, and
$E(B-V)_{\mathrm{max}}=0.5$. We finally apply attenutation through the IGM
following the prescriptions in Inoue et al. (2014). The observed $H_{160}$
magnitudes are pulled from a distribution constructed as the combination of a
uniform distribution in the range $22\leq H\leq 28$ and a truncated Gaussian
for fainter magnitudes with $\mu=26.5$, $\sigma=1$, $H_{\mathrm{min}}=22$, and
$H_{\mathrm{max}}=28$. We then normalize each spectral template to the
$H_{160}$ magnitude assigned to the simulated source.

Figure 15:  The fraction of real (grey histograms) and simulated (green
histograms) sources rejected in category 1e during visual inspection (Section
3.4). This category represents the only subjective rejections in our visual
classification scheme. We show the rejection fractions as a function of
$H_{160}$ magnitude (left panel), $H_{160}$ detection S/N (middle panel), and
half-light radius (right panel). The distributions of rejected real sources
approximately follow those for the rejected simulated sources, indicating a
bias likely introduced by our visual inspection that can directly affect our
measurement of the luminosity function. We use the left panel of this figure
to derive a correction factor that accounts for these rejection fractions as a
function of $H_{160}$ magnitude, and apply it to the calculated completeness
and effective volumes in Section 6.3.

We generate images of each source using Galfit (Peng et al., 2002) as Sérsic
profiles with index $n$ and axis ratio $b/a$ pulled randomly from truncated
Gaussian distributions with $\mu_{n}=1$ and $\sigma_{n}=2.5$ in the range
$0.5\leq n\leq 5$, and $\mu_{b/a}=0.8$ and $\sigma_{b/a}=0.2$, $0.1\leq
b/a\leq 1$, respectively. The half-light radii are determined by the absolute
magnitude of each source following the empirical size-luminosity relation
presented by Kawamata et al. (2018):

$r_{1/2}(M_{\mathrm{UV}})=r_{0}\times 10^{-0.4(M_{\mathrm{UV}}-M*)\beta},$ (3)

where $r_{0}$ is the modal radius at $M_{\mathrm{UV}}=-21$, $M^{*}$ is the
characteristic magnitude of the Schechter function fit to the luminosity
function, and $\beta$ is the slope of the size-luminosity relation. We adopt
the following parameter values from Table 2 of Kawamata et al. (2018):

$\displaystyle r_{0}/\mathrm{kpc}$ $\displaystyle=$
$\displaystyle\begin{cases}0.94&~{}~{}~{}~{}\text{for $z\leq 7.5$}\\\
0.58&~{}~{}~{}~{}\text{for $7.5<z\leq 8.5$}\\\ 0.42&~{}~{}~{}~{}\text{for
$z>8.5$}\end{cases}$ $\displaystyle\beta$ $\displaystyle=$
$\displaystyle\begin{cases}0.44&~{}~{}~{}~{}\text{for $z\leq 8.5$}\\\
0.40&~{}~{}~{}~{}\text{for $z>8.5$}\end{cases}$ (4) $\displaystyle M^{*}$
$\displaystyle=$ $\displaystyle\begin{cases}-20.56&\text{for $z\leq 7.5$}\\\
-19.95&\text{for $7.5<z\leq 8.5$}\\\ -19.80&\text{for $z>8.5$,}\end{cases}$

and additionally impose $\beta=0.25$ for sources with
$M_{\mathrm{UV}}\leq-21$, resulting in a size-luminosity relation described by
a broken power-law such that the slope flattens for brighter galaxies. In
order to account for intrinsic scatter in the relationship at fixed magnitude,
we apply a 0.2 dex scatter to the sizes derived from equation 3. We note that
the $r_{1/2}$ distribution measured in the Source Extractor catalogs of
recovered simulated sources traces that of the real sources, with the
simulated source sizes bracketing those of the real sources, confirming our
choice of parameter values for this relation. The Galfit model images for each
simulated source are convolved with the measured PSFs (Section 2.5)
corresponding to each dataset and filter.

Finally, we create 200 realizations of each of the 132 HST pointings
considered in this paper, each containing $\sim$100 sources randomly placed
across the image, for a total of $\sim$20000 input sources per pointing121212
As can be seen from Table 2, some fields cover significantly more or less area
than the average $\sim$4.6 arcmin2. In practice, we therefore add simulated
sources such that the density of sources in each pointing is constant at
$\sim$22 sources/arcmin2. This corresponds to 100 sources for each 4.6 arcmin2
field. However, we add fewer sources per realization and produce more
realizations for fields that have been partially masked as described in
Section 2.1. The larger fields are the result of combining two partially-
overlapping HST pointings into a single mosaic, and to these we add more
sources per realization. We perform at least 150 realizations per field, and
so the total number of simulated sources added to these larger fields exceeds
20000.. We then run our full sample selection process on these fields,
including source detection and photometry (Section 2.4), photometric redshift
fitting (Section 3.1), and the selection citeria used to select real high-
redshift candidates (Section 3.2). The one exception is that we do not add
Galactic extinction to the fluxes of the simulated sources, and so we skip the
Galactic extinction correction step described in Section 2.4. Finally, we
match the catalogs of recovered sources to the simulated input source
positions using a matching radius of 0$\farcs$5.

####  6.2.1 Quantifying Effects of Visual Inspection

As described in Section 3.4, the visual inspection process can affect sample
compeleteness by removing sources in a subjective way that is not replicable
with simulated sources. In an attempt to quantify this completeness effect, we
included the recovered simulated sources in our visual inspections, inspecting
an equal number of real and simulated sources. For each field, we randomly
selected the same number of simulated sources as there were real initial
candidates. In this way, the field-to-field filter coverages and noise
properties of the simulated sources matched those of the real candidates.
Importantly, this inspection was “anonymized”, or performed without knowing
which sources were simulated. This is the first time such a method – crucial
to correctly calculating sample completeness – has been used in high-redshift
galaxy selection.

Of the 166 real candidates that were rejected through visual inspection, we
focus here on the 81 rejections that were based on a subjective opinion
(category 1e, “a hot pixel, cosmic ray, or source with a strange morphology”,
Section 3.4). In Figure 15, we show the rejected fraction of the real and
simulated sources that were rejected in category 1e. The fraction of rejected
simulated sources increases towards fainter $H_{160}$ magnitude (left panel),
lower detection S/N in $H_{160}$ (middle), and smaller half-light radius
(right). Additionally, the fraction of rejected real sources largely follows
the same distribution as that for simulated sources, indicating that this
trend represents a bias in our visual inspection. However, we do not visually
inspect the many thousands of simulated sources that meet our selection
criteria. As we use the recovery of simulated sources to derive the effective
volume probed in each field, this bias imparted through our visual inspection
can therefore directly affect our measurement of the luminosity function. In
order to account for the rejected fraction of simulated sources, we use the
left panel of Figure 15 to apply a correction factor to the measured
completeness in bins of $H_{160}$ magnitude. For example, we decrease the
number of recovered sources in each field with $25.5<m_{H160}<26.0$ by 38%,
assuming that this percentage of sources would have been rejected during a
complete visual inspection. We discuss this correction further in Section 6.3.

Figure 16: The completeness and effective volume calculated for our selection
criteria. We plot the effective volume as a function of $H_{160}$ magnitude in
the left panel, showing each dataset separately with the total volume in
black. For all datasets, the lighter shaded region shows the effective volume
calculated directly from the completeness simulations, while the darker region
includes the correction for reduced completeness due to sources removed
through visual inspection. As expected given the distribution of rejected
fractions in the left panel of Figure 15, the effect of these corrections
increases towards fainter magnitudes. The thick black line indicates the
total, corrected volume for all three combined datasets. We show the results
of our completeness simulations as a function of redshift and $H_{160}$
magnitude in the right panel. The white circles indicate the redshifts and
magnitudes of our 13 candidates. The completeness is highest for bright
sources and falls off towards fainter magnitude. It also peaks at $z>10$, when
the Ly$\alpha$ break redshifts out of the $J$ band and sources are
$H_{160}$-only single detections. At these redshifts, non-detections in all
three blue filters result in larger portions of the redshift PDF at $z>8$, and
the recovery fraction in the WISP fields increases significantly and boosts
the overall completeness.

###  6.3 Effective Volume

We calculate the number density of UV-bright galaxies at $z\sim 9$ using the
effective volume method (Felten, 1977). This method involves determining the
maximum volume within which a galaxy of a given luminosity would be detected
by these HST observations and selected with our specific set of criteria. The
effective volume, $V_{\mathrm{eff}}$, is given by

$V_{\mathrm{eff}}(M_{\mathrm{UV}})=\int\frac{\mathrm{d}V}{\mathrm{d}z}P(M_{\mathrm{UV}},z)\>\mathrm{d}z,$
(5)

where $\mathrm{d}V/\mathrm{d}z$ is the differential comoving volume element,
and $P(M_{\mathrm{UV}},z)$ is the probability that a source of absolute
magnitude $M_{\mathrm{UV}}$ and redshift $z$ will make it into our sample
given our selection function. We estimate $P(M,z)$ using the set of field-
specific simulations described in Section 6.2.

We calculate the completeness as the ratio of the number of recovered sources
to input sources in bins of input magnitude and redshift. In order to account
for the simulated sources that we rejected during our visual inspection, we
decrease the number of recovered sources in each magnitude bin by the fraction
of rejected simulated sources shown in Figure 15. For this correction, we
scale the magnitude bins to the 5$\sigma$ $H_{160}$ magnitude limit of each
field. In this way we approximate rejecting an increasing number of simulated
sources as we approach the field depth. We do not include the non-zero
fraction of bright ($m_{H160}\sim 22.3$) rejected simulated sources, as there
are no rejected real sources in that magnitude bin.

We show the resulting completeness as a function of redshift and $H_{160}$
magnitude in the right panel of Figure 16. There are almost no recovered
sources at $z<7.5$, where the Ly$\alpha$ break falls close to the center of
the $Y_{098}$ and $Y_{105}$. In these cases, the sources are detected in these
filters with S/N$>$2 and therefore do not satisfy our selection criteria. As
the $Y_{098}$ is narrower than the $Y_{105}$, simulated sources at $z\sim 7.5$
are recovered in BoRG[$z$8] fields, while the recovery fraction in HIPPIES GO
12286 fields starts to increase for $z>8$. The completeness dips slightly at
$z\sim 8.6$, the redshift at which the Ly$\alpha$ break transitions almost
entirely out of the $Y_{105}$ and sources in BoRG[$z$8] and HIPPIES fields
become two-band detections. There is another dip at $z\sim 10$, as the
Ly$\alpha$ break shifts out of the $J_{125}$ and $J_{110}$ filters and sources
are detected in $H_{160}$ only. After each initial dip, the completeness
increases with redshift as more filters lie blueward of the break, ruling out
lower-redshift solutions and increasing the portion of the redshift PDF that
lies at $z>8$. The recovery in WISP fields is very low at $z\lesssim 10$
because the broad $J_{110}$ filter results in poorly constrained photometric
redshift PDFs for sources from $z\sim 7$ to $z\sim 10$. The completeness is
therefore highest at $z>10.5$, when the WISP recovery fractions increase again
for $H_{160}$-only detections.

We calculate the effective volume probed in each field using equation 5. In
the left panel of Figure 16, we show the effective volume calculated as a
function of $H_{160}$ apparent magnitude. For each dataset, the darker shaded
region shows the $V_{\mathrm{eff}}$ that has been corrected for the bias
introduced by our visual inspections, and the lighter shaded region shows the
original, uncorrected $V_{\mathrm{eff}}$. It can be seen here how this
correction increasingly reduces the effective volume for fainter $H_{160}$
magnitudes. The thick solid black line indicates the total, corrected volume
of all three combined datasets. We use the effective volume calculated as a
function of absolute UV magnitude for our calculation of number densities in
Section 6.4, adopting a rest-wavelength of 1500Å that corresponds to
approximately the center of the $H_{160}$ filter at $z\sim 9.5$.

###  6.4 Volume Number Densities

We calculate the observed number density of galaxies in the redshift range
$8.5\lesssim z\lesssim 11$ following the methods described in Finkelstein et
al. (2022). We describe our steps here but refer the reader to Finkelstein et
al. (2022) for more details. We use a Markov Chain Monte Carlo (MCMC) analysis
with the Cash statistic (Cash, 1979; Ryan et al., 2011) as the goodness-of-fit
measure. For a given assumed number density, the Cash statistic describes the
likelihood that the observed number of galaxies matches the expected number,
but models the probability distribution as a Poisson rather than a Gaussian
distribution. This choice of statistic is therefore appropriate for our small
sample size.

Additionally, we use the “pseudo-binning” technique developed by Finkelstein
et al. (2022), which reduces the biases imparted by choice of absolute
magnitude bin center and width on the recovered shape of the luminosity
function. At each step of the MCMC chain, the number density is calculated for
a given $M_{\mathrm{UV}}$ in a bin with a width randomly pulled from the range
$0.2-1.5$ mag. The volume associated with this bin is determined using the
completeness simulations described in Section 6.3 and equation 5. We pull an
absolute magnitude for each source randomly from the distributions calculated
in Section 6.1, using the magnification-corrected distributions for the three
strongly-lensed sources. We then determine the observed number of sources in
the bin given these randomly-drawn quantities. We perform this calculation in
magnitude steps of $\Delta M_{\mathrm{UV}}=0.1$ in the range
$-24<M_{\mathrm{UV}}<-20$. In this way, the number density posterior
distribution incorporates the uncertainties on the absolute magnitude of each
candidate – which in turn fold in the uncertainties in redshift, $H_{160}$
observed magnitude, and magnification correction – as well as reducing the
dependence on choice of magnitude bins. We use 20 walkers with a burn in of
$10^{5}$ steps, after which we discard the chains and use an additional
$10^{4}$ steps to calculate the posterior distribution.

We present the calculated number densities in Table 7 in the Appendix and
display our results in Figure 17, plotting the 68% and 95% ranges of the
density posterior as blue shaded regions. The shading extends from
$M_{\mathrm{UV}}=-23.05$, the brightest bin that contains a candidate, to
$M_{\mathrm{UV}}=-20.5$, the absolute magnitude at which our average
completeness drops below 20%. We plot vertical blue error bars throughout the
shaded region to indicate the uncertainty in our measurements due to cosmic
variance. We estimate this uncertainty using the calculator from Bhowmick et
al. (2020). At all magnitudes, the uncertainty on the number density due to
cosmic variance is very low. The overall uncertainties are dominated by
Poisson noise. For comparison, we also show a number of results from the
literature, including the pseudo-binning results from Finkelstein et al.
(2022) for their 11 candidates in the same redshift range (purple shaded
region). We discuss the implications of our results in the context of other
similar high-redshift searches in Section 7.

Figure 17: The UV luminosity function (blue shaded region) in the range
$8.5\lesssim z\lesssim 11$ as measured in the BoRG[$z$8], HIPPIES GO 12286 and
WISP datasets using the pseudo-binning technique developed by Finkelstein et
al. (2022). The inner (outer) blue shaded region is the Poisson uncertainty
determined by the 68% (95%) range of the number density posterior, extending
from the brightest bin that contains a candidates to the magnitude at which
our average completeness drops below 20%. The solid blue line shows the upper
68% range at brighter and fainter magnitudes, and the blue vertical errorbars
indicate the uncertainty due to cosmic variance based on the calculator from
Bhowmick et al. (2020). We show the absolute magnitudes of each of our 13
candidates at the top as blue circles, with errorbars representing their
$M_{\mathrm{UV}}$ posterior distribution 68% ranges. We compare our results
with those from Finkelstein et al. (2022) (purple shaded region), as well as
other studies using HST parallel fields (Bernard et al., 2016; Calvi et al.,
2016; Livermore et al., 2018; Morishita et al., 2018; Rojas-Ruiz et al.,
2020), ground-based imaging in the COSMOS/UltraVISTA field (Stefanon et al.,
2019; Bowler et al., 2020), and HST imaging in the extragalactic legacy fields
(Oesch et al., 2018; Bouwens et al., 2021). For reference, we also show the
Schechter function fits from Oesch et al. (2018) and Bouwens et al. (2021),
the Schechter function predictions from Finkelstein (2016), and the double
power law luminosity functions from Bowler et al. (2020). We find a high
density of bright galaxies, in good agreement with other HST parallel studies.
Our luminosity function is consistent for $M_{\mathrm{UV}}\lesssim-21.5$ with
that of Finkelstein et al. (2022), which was calculated using the same method
but in a handful of fields. The uncorrelated nature of the pure-parallel
observations we include in our analysis significantly reduces the effects of
cosmic variance, allowing for a more complete sampling of the population of
bright galaxies at $z\sim 9-10$.

##  7 Discussion

###  7.1 Comparison with Previous Luminosity Function Calculations

We begin by comparing our results to other studies in similar redshift ranges.
In Figure 17, we include number density measurements from high-redshift galaxy
searches performed in HST parallel fields (Bernard et al., 2016; Calvi et al.,
2016; Livermore et al., 2018; Morishita et al., 2018; Rojas-Ruiz et al.,
2020), ground-based imaging in the COSMOS/UltraVISTA field (Stefanon et al.,
2019; Bowler et al., 2020), the CANDELS fields (Finkelstein et al., 2022) as
well as those that included the Hubble Frontier Fields (Oesch et al., 2018)
and the Hubble Ultra Deep Field (Bouwens et al., 2021). Collectively, these
surveys cover a large dynamic range in both magnitude and volume, yet
individually they have complementary strengths. The Hubble extragalactic
legacy fields have the depth needed to probe to fainter magnitudes, yet
include only a handful of lines of sight and are therefore more susceptible to
cosmic variance. On the other hand, the HST pure parallel programs include
hundreds of independent pointings, and the ground-based COSMOS/UltraVISTA
field covers a very large area at a shallow depth. These surveys are more
likely to detect the rare, brightest galaxies, but at the expense of lower
sample completeness for $M_{\mathrm{UV}}\gtrsim-21.5$.

Our results are consistent with other searches performed in HST parallel
surveys. Calvi et al. (2016), Morishita et al. (2018) and Rojas-Ruiz et al.
(2020) measured the rest-UV luminosity function in the BoRG[$z910$] survey,
which includes imaging in $JH_{140}$. The second detection filter helps
protect against the selection of spurious sources and improves the photometric
redshift by more fully constraining the location of the Ly$\alpha$ break. The
132 HST fields we include in our results are less optimized for the selection
of $z\sim 9-10$ galaxies. Nevertheless, the number densities from these works
at $z\sim 9$ (grey symbols), $z\sim 10$ (orange symbols), and across the full
redshift range $8.4<z<10.6$ as in Rojas-Ruiz et al. (2020, orange stars) lie
approximately at the median density we calculate across the full absolute
magnitude range that we consider. Bernard et al. (2016) explore some of the
same fields that we do, yet our two samples of high-redshift galaxy candidates
do not overlap. However, the number densities reported by Bernard et al.
(2016) at $z\sim 10$ are in good agreement with our measurements. We note that
while the UV luminosity function from the SuperBoRG fields has not been
published at the time of this writing, Roberts-Borsani et al. (2021) report
the identification of 49 candidates in the range $z\sim 8-12$, likely
resulting in a high number density in agreement with these other searches in
pure-parallel data.

Overall, the results using the HST pure parallel surveys find number densities
above those measured in the HST legacy or ground-based imaging. The two
brightest magnitude bins at $z\sim 9$ from Bouwens et al. (2021) are
consistent with our measurement within the uncertainties, yet the Oesch et al.
(2018) measurements at $z\sim 10$ are much lower. Both studies include fields
deep enough to probe the $z\sim 9-10$ galaxy population down to
$M_{\mathrm{UV}}>-18$, yet do not find galaxies at $M_{\mathrm{UV}}<-22$. The
brightest $z\sim 10$ candidate contributing to the luminosity function of
Oesch et al. (2018) is GN-z11 (Oesch et al., 2016) with a magnitude of
$m_{\mathrm{UV}}=-21.6$. The bright end of the luminosity function measured by
Stefanon et al. (2019) is also lower but consistent with our measurements
within the uncertainties.

We next focus on the comparison between our results and those of Finkelstein
et al. (2022), as we have performed very similar sample selections,
completeness analyses, and have used the same method to measure the UV
luminosity function. Our luminosity functions are fully consistent down to
$M_{\mathrm{UV}}\sim-21.8$, below which the purple region decreases and is
more consistent with the Schechter functions of Finkelstein (2016, predicted
for a smooth decline in the luminosity function at $z>8$), Bowler et al.
(2020, showing a slower evolution in the bright end of the luminosity function
at $z>8$) and Bouwens et al. (2021, proposing an accelerated decline in the
luminosity function). The high-redshift sample presented by Finkelstein et al.
(2022) include an overdensity in the EGS field, which may serve to increase
their number densities over similar studies in the HST legacy fields. At the
bright end, however, the uncertainties of Finkelstein et al. (2022) are
dominated by cosmic variance, finding a fractional uncertainty of 0.95 at
$M_{\mathrm{UV}}=-22.5$. With over 130 independent pointings included in our
analysis, we find a fractional uncertainty of 0.14 at $M_{\mathrm{UV}}=-22.8$
(blue vertical errorbars). This indicates that the high density of bright
galaxies as reported here and in Finkelstein et al. (2022) is robust against
cosmic variance, though we note that the Poisson uncertainties in our
measurements are large. As described in Section 5.2, we expect our
contamination from lower-redshift interlopers to increase towards fainter
magnitudes. This increasing contamination could explain why we calculate a
higher luminosity function than that presented by Finkelstein et al. (2022)
for $M_{\mathrm{UV}}\gtrsim-21.8$.

These results indicate that wide-area or pure parallel programs provide
valuable insight about the bright end of the rest-UV luminosity function.
Their large areas or many uncorrelated pointings significantly reduce the
effects of cosmic variance and allow for a more complete sampling of the
population of bright, $z\sim 9-10$ galaxies. While there is significant
scatter and Poisson uncertainty for $M_{\mathrm{UV}}\lesssim-22$, there is a
broad consensus among these programs that there is an excess of bright
galaxies at $z\sim 9-10$ compared to a Schechter function parameterization of
the luminosity function. By $M_{\mathrm{UV}}\sim-21.5$, results from these
shallow observations may be more prone to contamination, and we should turn to
deeper fields for measurements of the luminosity function.

While many factors affect the evolving abundance of galaxies with redshift,
the two dominant processes influencing the bright end of the observed UV
luminosity function at $z\sim 9$ are likely star formation efficiencies and
dust attenuation. Somerville et al. (2015) and Yung et al. (2019a) showed that
the predicted number density of high redshift galaxies ($z\gtrsim 4$) is very
sensitive to the efficiency with which molecular hydrogen is converted into
stars. This is because at these redshifts, the molecular gas consumption time
(a few Gyr in nearby galaxies) becomes comparable to the age of the Universe.
Specifically, in their fiducial model, Yung et al. (2019a) adopt a double
power law for the relationship between molecular gas surface density and star
formation rate surface density, where the slope of the power law becomes
steeper above a critical molecular gas surface density. Galaxies at early
times are more compact than those at lower redshifts, and so most star
formation is likely occurring in gas at densities above this critical density.
The net effect is that gas depletion times are shorter, and star formation
efficiencies higher, at high redshift. The authors find that with a slope of
$\sim 2$ at high gas densities (a bit steeper than the slope of unity measured
for nearby spiral galaxies; Bigiel et al., 2008), their semi-analytic model
can reproduce the observed number density of bright galaxies out to $z\sim 8$.
A model with a single power law star formation relation with a slope of either
unity or 1.5 (as in the canonical Kennicutt-Schmidt relation Schmidt, 1959,
1963; Kennicutt, 1989, 1998) significantly underproduces bright galaxies even
by $z\sim 6$.

Dust attenuation can also play a major role in determining the shape of the
bright end of the UV luminosity function (e.g., Vogelsberger et al., 2020).
The high density of UV-bright galaxies in Figure 17 may indicate low or
negligible dust content in massive galaxies at these redshifts. The shallow
evolution in the bright end at high redshift may be cause by a redshift and/or
mass dependence in dust production and dust-to-metal ratios. Finkelstein et
al. (2022) shows that the semi-analytic model from Yung et al. (2019a) can
best reproduce their observed density of bright galaxies at $z\sim 9$ when
there is no dust attenuation included (yet SED fitting for the galaxies in
their sample indicate non-negligible dust contents Tacchella et al., 2022). It
is also important to note that the effects on the bright end of lower dust
attenuation are degenerate with those of higher star formation efficiencies,
and detailed observations with JWST NIRCam and MIRI will be required to
separate these effects in bright, high-redshift galaxies.

This observed excess could also be due to a high fraction of faint active
galactic nuclei (AGNs) in our sample. While there have been only a handful of
AGN detected at $z>7$, massive accreting black holes must exist at earlier
times to explain the masses of the quasars observed at $z\sim 6-7$ (Bañados et
al., 2018). Additionally, there is a growing body of evidence that some of the
brightest known $z>7$ star-forming galaxies may in fact host detectable AGN
activity, typically detected via significant N V emission (e.g., Tilvi et al.,
2016; Laporte et al., 2017; Hu et al., 2017) including one of the highest-
redshift spectroscopically-confirmed galaxies at $z=8.68$ (Zitrin et al.,
2015; Mainali et al., 2018). We cannot rule out the possibility that the
luminosity of these bright galaxies is partially due to an AGN contribution.
Spectroscopic followup with NIRSpec’s medium and high resolution gratings
could help shed light on this question.

###  7.2 Bright Galaxies during the Epoch of Reionization with JWST

This sample of 13 galaxy candidates at $8.5\lesssim z\lesssim 11$ represent
the very bright end of the rest-UV luminosity function, with
$M_{\mathrm{UV}}\lesssim-21.2$. Such massive galaxies likely formed in regions
of relative overdensity, where the build-up of mass began at earlier times.
These regions are of particular interest in understanding the early stages of
the reionization of the IGM. Whether the progression of the epoch of
reionization at early times was dominated by ionizing emission from massive
galaxies with high escape fractions (e.g., Naidu et al., 2020, 2021; Matthee
et al., 2021), or from the more numerous lower-mass population (e.g.,
Finkelstein et al., 2019), these overdense regions would have been among the
first to ionize. Bright galaxies may therefore reside in ionized bubbles even
as early as $z\gtrsim 8$. Indeed, Ly$\alpha$ has been detected in the spectra
of two galaxies at $z\sim 8.7$ in the EGS field (Zitrin et al., 2015; Larson
et al., 2022), potentially indicating ionized regions large enough for the
Ly$\alpha$ photons to redshift out of resonance before encountering the
neutral IGM.

The $z\sim 9-10$ candidates we present in this paper may trace overdense
regions at early times, making them excellent targets for follow-up with JWST.
Deep NIRCam imaging can be used to detect the population of fainter galaxies
as demonstrated by Castellano et al. (2016, 2018, at $z\sim 7$), and deep
NIRSpec spectroscopy can be used to probe reionization in these regions via
Ly$\alpha$ and Mg II (e.g., Henry et al., 2018). The NIRSpec prism also
provides the opportunity to characterize the ionizing efficiencies of bright
galaxies. The building narrative surrounding reionization hints at a chaotic
early universe, with galaxies undergoing intense bursts of star-formation,
producing massive, bright stars that power strong nebular emission lines
(e.g., Stark, 2016; Finkelstein, 2016; Endsley et al., 2020). The relative
strength and ratios of these lines, such as that of [O III]$\lambda\lambda
4959,5007$/[O II]$\lambda\lambda 3727,3729$ (=O32; e.g., Steidel et al.,
2016), can therefore reveal significant detail about the strength of the
ionizing continuum in galaxies at $z>9$ (e.g., Shapley et al., 2003; Erb et
al., 2010; Steidel et al., 2016; Ravindranath et al., 2020). Observational
evidence shows that high O32 correlates with an increasing ionizing efficiency
(e.g., Tang et al., 2018). For galaxies during the epoch of reionization,
large measured ionizing efficiencies may indicate sources that are capable of
ionizing a big bubble in the IGM. The combination of NIRSpec spectroscopy and
deep NIRCam imaging will test whether reionization is driven by the most
massive galaxies, or by the smaller, fainter sources.

With JWST GO-2426 (Co-PIs: Bagley & Rojas-Ruiz), we will follow-up five
sources from this paper (Par0456-2203_473, Par0756+3043_92, Par0956-0450_684,
Par335_251, and Par2346-0021_164131313While we removed Par2346-0021_164 from
consideration when calculating the luminosity function at $z\sim 8.5-10$,
$>$60% of the integrated redshift PDF is at $z>8$. This candidate may
therefore still be a high-redshift galaxy. Alternatively, spectroscopic
characterization of a low-redshift interloper is also valuable in evaluating
high-redshift galaxy selection techniques.) as well as six additional
candidates presented in Rojas-Ruiz et al. (2020). With approximately one hour
of exposure time per target with the NIRSpec prism in fixed slit mode, the
observations will be sensitive enough to detect the Ly$\alpha$ break for
redshift confirmation as well as nebular emission lines such as [O II],
H$\beta$, and [O III]. There are a number of other JWST programs that will be
exploring bright galaxies during the epoch of reionization, including GO 1740
(PI: Harikane), GO 1747 (PI: Roberts-Borsani), GO 1933 (PI: Matthee), GO 2279
(PI: Naidu), and GO 2659 (PI: Weaver). Together, these Cycle 1 JWST
observations will significantly expand our understanding of massive galaxies
at early times and their contribution to the reionization of the IGM.

##  8 Summary

We have presented a search for bright $z\sim 9-10$ candidate galaxies in 132
HST fields covering 620 arcmin2. The fields were observed in parallel as part
of the BoRG[$z$8], HIPPIES GO 12286, and WISP surveys. The independent and
uncorrelated nature of these pointings significantly reduces the effect of
cosmic variance compared to observations covering larger contiguous areas.
While the three datasets we include in our analysis use slightly different
filter sets, all provide imaging in $H_{160}$, a $J$-band filter (either
$J_{125}$ or $J_{110}$), and a subset of $V_{606}$, $I_{600}$, $I_{814}$,
$Y_{098}$ and $Y_{105}$. This filter suite allows for two-band detections of
galaxies at $z\sim 9-10$ and $H_{160}$-only detections at $z\gtrsim 10$. The
majority of the HST pointings were also observed by Spitzer/IRAC at 3.6µm and
(for a smaller subset) 4.5µm. The HST imaging in two bluer filters is needed
to detect the redshifted Ly$\alpha$ break, and the IRAC imaging at longer
wavelengths is crucial for ruling out contamination by passive or dusty star-
forming galaxies at $z\sim 2-3$ and low mass stars or brown dwarfs.

We performed a careful $z\gtrsim 8.5$ sample selection using criteria
involving detection significance (S/N${}_{160}>5$) and photometric redshift
fitting. We used EAZY to measure photometric redshifts and selected sources
with at least 70% of their redshift probability distribution function at
$z>8$. We also enforced a lower limit for source half-light radius, using a
random forest algorithm to set the radius thresholds for each dataset. We
conducted a thorough check for image persistence, exploring all observations
obtained in the 24 hours before each image of a candidate field. We performed
an additional screening for self-persistence in WISP images that is caused by
zeroth orders in the grism exposures.

Even the most conservative selection criteria cannot remove all spurious
detections. We therefore visually inspected all candidates, rejecting
diffraction spikes, sources along image edges, bad pixels, and sources with
morphologies that appeared similar to hot pixels, cosmic rays, or other
artifacts. Importantly, we performed this visual inspection on both real and
simulated candidates, randomizing their order and anonymizing their nature.
This technique allowed us to include the effects of these subjective
rejections in our completeness corrections.

We next measured the photometry of each candidate in all available
Spitzer/IRAC imaging, using Galfit to model the light from the candidates and
any neighboring sources in the IRAC images. We refit the HST+Spitzer
photometry with EAZY and rejected one source with a redshift probability
distribution that no longer satisfied our selection criteria. Additionally, we
rejected one candidate with a low $\Delta\chi^{2}$ that indicated the high-
redshift solution was not preferred at high significance. Finally, we explored
whether any of the candidates were likely to be stars by fitting the
HST+Spitzer photometry with spectral models of M, L, and T dwarf stars. All
candidates were better fit by a high-redshift galaxy template. Our final
sample was thirteen candidates from $8.3\lesssim z\lesssim 11$, including the
first $z>8$ candidate identified from the WISP survey.

We presented our sample in Section 3, including a comparison to previous
studies in the same fields. While our selection criteria are not sensitive to
the $z\sim 7-8$ samples previously identified in the BoRG[$z$8] fields, we did
recover two candidates presented by Bradley et al. (2012) that likely lie at
the high end of their redshift distribution. We did not recover any of the
three $z\sim 9-10$ candidates presented by Bernard et al. (2016), highlighting
the dependence of high-redshift candidate samples on the methods used to
select them. We then explored the possibility of low-redshift contamination in
our high redshift sample in two ways. First, we stacked the imaging in all
filters and did not detect significant flux in any of the dropout bands. Next,
we estimated the contamination from faint low-redshift galaxies by dimming the
fluxes of bright $z=1-4$ galaxies and reselecting them with our selection
criteria. We found low contamination rates of $\sim$1-3% in almost all fields
(with the highest rate of 12.9% in the WISP field).

After applying a correction for magnification along the line of sight to the
three sources with bright, close neighbors, we calculated the UV absolute
magnitudes of each candidate. Our magnitude uncertainties included the
uncertainties in the photometric redshifts, $H_{160}$ fluxes, and
magnification corrections. We used simulated sources to quantify our sample
selection incompleteness and determine the effective volume of the combined
132 pointings. We then calculated the rest-UV luminosity function using an
MCMC analysis that incorporates our magnitude uncertainties and the pseudo-
binning technique presented by Finkelstein et al. (2022). We found number
densities that are higher than those presented in most searches in HST legacy
fields, but are in good agreement with other searches performed in HST
parallel fields. Our luminosity function is also consistent with that of
Finkelstein et al. (2022) down to $M_{\mathrm{UV}}\sim-21.8$. Taken together,
these results demonstrate the power of pure parallel programs in measuring the
bright end of the luminosity function with a much smaller uncertainty due to
cosmic variance. However, for fainter magnitudes, the shallow parallel
observations may be more prone to contamination.

These results also indicate evidence for an excess of bright
($M_{\mathrm{UV}}\sim-22$) $z\sim 9-10$ galaxies compared to expectations for
a smoothly evolving Schechter function from lower redshifts. We concluded by
discussing possible implications, including the potential for increased star
formation efficiencies and decreased dust attenutation to explain the
observations. We also discussed what can be learned by studying bright
galaxies during the epoch of reionization. These galaxies likely trace
overdense regions that may have been among the first to ionize. They are also
excellent targets for efficient followup with JWST, bright enough for detailed
spectroscopic analysis of their ionizing power and contribution to
reionization. Five of the sources presented in this paper will be observed
with the NIRSpec prism as part of JWST GO-2426, and many other bright high-
redshift galaxies will be observed during Cycle 1. Our understanding of the
early universe stands to be transformed in the coming year of observations.

MBB, SLF, and KDF acknowledge that they work at an institution, the University
of Texas at Austin, that sits on indigenous land. The Tonkawa lived in central
Texas and the Comanche and Apache moved through this area. We pay our respects
to all the American Indian and Indigenous Peoples and communities who have
been or have become a part of these lands and territories in Texas. We are
grateful to be able to live, work, collaborate, and learn on this piece of
Turtle Island.

In this paper we refer to the James Webb Space Telescope using only the
acronym JWST, reflecting our choice to celebrate the promise of this telescope
without acknowledging the public official for whom it is named. This
individual has been implicated in anti-LGBTQI+ attitudes that do not reflect
the authors’ values related to inclusion in science.

The authors are grateful to Mira Mechtley for developing the comprehensive
pipeline use to reduce the HIPPIES GO 12286 dataset presented in this paper.

This research was partially supported by NASA through ADAP awards NNX16AN47G
(linked to Spitzer GO program 11121) and 80NSSC18K0954. MBB acknowledges
support from a NASA Keck PI Data Award, PID 83/2019B_N127, administered by the
NASA Exoplanet Science Institute. SRR acknowledges financial support from the
International Max Planck Research School for Astronomy and Cosmic Physics at
the University of Heidelberg (IMPRS–HD). RSS is supported by the Simons
Foundation. YSD acknowledges the support from the NSFC grants 11933003, and
the China Manned Space Project with No. CMS-CSST-2021-A05.

This research is based on observations made with the NASA/ESA Hubble Space
Telescope obtained from the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
under NASA contract NAS 5–26555. These observations are associated with
programs 11519, 11520, 11524, 11528, 11530, 11533, 11534, 11541, 11700, 11702,
12024, 12025, 12283, 12286, 12572, 12902, 12905, 13352, 13517 and 14178.

This research has made use of the NASA/IPAC Infrared Science Archive, which is
funded by the National Aeronautics and Space Administration and operated by
the California Institute of Technology.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Atek et al. (2010) Atek, H., Malkan, M., McCarthy, P., et al. 2010, ApJ, 723, 104, doi: 10.1088/0004-637X/723/1/104
  * Atek et al. (2011) Atek, H., Siana, B., Scarlata, C., et al. 2011, ApJ, 743, 121, doi: 10.1088/0004-637X/743/2/121
  * Atek et al. (2015) Atek, H., Richard, J., Kneib, J.-P., et al. 2015, ApJ, 800, 18, doi: 10.1088/0004-637X/800/1/18
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * Beckwith et al. (2006) Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729, doi: 10.1086/507302
  * Beers et al. (1990) Beers, T. C., Flynn, K., & Gebhardt, K. 1990, AJ, 100, 32, doi: 10.1086/115487
  * Behroozi et al. (2020) Behroozi, P., Conroy, C., Wechsler, R. H., et al. 2020, arXiv e-prints, arXiv:2007.04988.  https://arxiv.org/abs/2007.04988
  * Behroozi & Silk (2015) Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32, doi: 10.1088/0004-637X/799/1/32
  * Bernard et al. (2016) Bernard, S. R., Carrasco, D., Trenti, M., et al. 2016, ApJ, 827, 76, doi: 10.3847/0004-637X/827/1/76
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astronomy and Astrophysics Supplement Series, 117, 393, doi: 10.1051/aas:1996164
  * Bhowmick et al. (2020) Bhowmick, A. K., Somerville, R. S., Di Matteo, T., et al. 2020, MNRAS, 496, 754, doi: 10.1093/mnras/staa1605
  * Bigiel et al. (2008) Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846, doi: 10.1088/0004-6256/136/6/2846
  * Blanc et al. (2008) Blanc, G. A., Lira, P., Barrientos, L. F., et al. 2008, ApJ, 681, 1099, doi: 10.1086/588018
  * Bouwens et al. (2017) Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017, ApJ, 843, 129, doi: 10.3847/1538-4357/aa70a4
  * Bouwens et al. (2019) Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25, doi: 10.3847/1538-4357/ab24c5
  * Bouwens et al. (2009) Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936, doi: 10.1088/0004-637X/705/1/936
  * Bouwens et al. (2010) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2010, ApJ, 709, L133, doi: 10.1088/2041-8205/709/2/L133
  * Bouwens et al. (2015) —. 2015, ApJ, 803, 34, doi: 10.1088/0004-637X/803/1/34
  * Bouwens et al. (2016) Bouwens, R. J., Oesch, P. A., Labbé, I., et al. 2016, ApJ, 830, 67, doi: 10.3847/0004-637X/830/2/67
  * Bouwens et al. (2021) Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, arXiv e-prints, arXiv:2102.07775.  https://arxiv.org/abs/2102.07775
  * Bower et al. (2012) Bower, R. G., Benson, A. J., & Crain, R. A. 2012, MNRAS, 422, 2816, doi: 10.1111/j.1365-2966.2012.20516.x
  * Bowler et al. (2020) Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059
  * Bowler et al. (2014) Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2014, MNRAS, 440, 2810, doi: 10.1093/mnras/stu449
  * Bowler et al. (2015) —. 2015, MNRAS, 452, 1817, doi: 10.1093/mnras/stv1403
  * Bradley et al. (2019) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2019, astropy/photutils: 1.0.0, 0.6, Zenodo, doi: 10.5281/zenodo.2533376.  https://doi.org/10.5281/zenodo.2533376
  * Bradley et al. (2012) Bradley, L. D., Trenti, M., Oesch, P. A., et al. 2012, ApJ, 760, 108, doi: 10.1088/0004-637X/760/2/108
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Bridge et al. (2019) Bridge, J. S., Holwerda, B. W., Stefanon, M., et al. 2019, ApJ, 882, 42
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000, doi: 10.1046/j.1365-8711.2003.06897.x
  * Burgasser (2014) Burgasser, A. J. 2014, in Astronomical Society of India Conference Series, Vol. 11, Astronomical Society of India Conference Series, 7–16
  * Calvi et al. (2016) Calvi, V., Trenti, M., Stiavelli, M., et al. 2016, ApJ, 817, 120, doi: 10.3847/0004-637X/817/2/120
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Casertano et al. (2000) Casertano, S., de Mello, D., Dickinson, M., et al. 2000, AJ, 120, 2747, doi: 10.1086/316851
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Castellano et al. (2016) Castellano, M., Dayal, P., Pentericci, L., et al. 2016, ApJ, 818, L3, doi: 10.3847/2041-8205/818/1/L3
  * Castellano et al. (2018) Castellano, M., Pentericci, L., Vanzella, E., et al. 2018, ApJ, 863, L3, doi: 10.3847/2041-8213/aad59b
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Clampin et al. (2000) Clampin, M., Ford, H. C., Bartko, F., et al. 2000, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4013, Advanced camera for surveys, ed. J. B. Breckinridge & P. Jakobsen, 344–351
  * Coe et al. (2015) Coe, D., Bradley, L., & Zitrin, A. 2015, ApJ, 800, 84, doi: 10.1088/0004-637X/800/2/84
  * Coe et al. (2013) Coe, D., Zitrin, A., Carrasco, M., et al. 2013, ApJ, 762, 32, doi: 10.1088/0004-637X/762/1/32
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Ellis et al. (2001) Ellis, R., Santos, M. R., Kneib, J.-P., & Kuijken, K. 2001, ApJ, 560, L119, doi: 10.1086/324423
  * Endsley et al. (2020) Endsley, R., Stark, D. P., Chevallard, J., & Charlot, S. 2020, arXiv e-prints, arXiv:2005.02402.  https://arxiv.org/abs/2005.02402
  * Erb et al. (2010) Erb, D. K., Pettini, M., Shapley, A. E., et al. 2010, ApJ, 719, 1168, doi: 10.1088/0004-637X/719/2/1168
  * Faisst et al. (2016) Faisst, A. L., Capak, P., Hsieh, B. C., et al. 2016, ApJ, 821, 122, doi: 10.3847/0004-637X/821/2/122
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10, doi: 10.1086/422843
  * Felten (1977) Felten, J. E. 1977, AJ, 82, 861, doi: 10.1086/112140
  * Finkelstein et al. (2015a) Finkelstein, K. D., Finkelstein, S. L., Tilvi, V., et al. 2015a, ApJ, 813, 78, doi: 10.1088/0004-637X/813/1/78
  * Finkelstein (2016) Finkelstein, S. L. 2016, Publications of the Astronomical Society of Australia, 33, e037, doi: 10.1017/pasa.2016.26
  * Finkelstein et al. (2010) Finkelstein, S. L., Papovich, C., Giavalisco, M., et al. 2010, ApJ, 719, 1250, doi: 10.1088/0004-637X/719/2/1250
  * Finkelstein et al. (2012a) Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012a, ApJ, 758, 93, doi: 10.1088/0004-637X/758/2/93
  * Finkelstein et al. (2012b) Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012b, ApJ, 756, 164, doi: 10.1088/0004-637X/756/2/164
  * Finkelstein et al. (2015b) Finkelstein, S. L., Ryan, Russell E., J., Papovich, C., et al. 2015b, ApJ, 810, 71, doi: 10.1088/0004-637X/810/1/71
  * Finkelstein et al. (2015c) Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015c, ApJ, 814, 95, doi: 10.1088/0004-637X/814/2/95
  * Finkelstein et al. (2019) Finkelstein, S. L., D’Aloisio, A., Paardekooper, J.-P., et al. 2019, ApJ, 879, 36, doi: 10.3847/1538-4357/ab1ea8
  * Finkelstein et al. (2022) Finkelstein, S. L., Bagley, M., Song, M., et al. 2022, ApJ, 928, 52, doi: 10.3847/1538-4357/ac3aed
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63, doi: 10.1086/316293
  * Fowler et al. (2017) Fowler, J., Markwardt, L., Bourque, M., & Anderson, J. 2017, Analysis of Dragon’s Breath and Scattered Light Detector Anomalies on WFC3/UVIS, Space Telescope WFC Instrument Science Report
  * Froning & Green (2009) Froning, C. S., & Green, J. C. 2009, Ap&SS, 320, 181, doi: 10.1007/s10509-008-9758-y
  * Gnedin (2016) Gnedin, N. Y. 2016, ApJ, 825, L17, doi: 10.3847/2041-8205/825/2/L17
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35, doi: 10.1088/0067-0049/197/2/35
  * Henry et al. (2018) Henry, A., Berg, D. A., Scarlata, C., Verhamme, A., & Erb, D. 2018, ApJ, 855, 96, doi: 10.3847/1538-4357/aab099
  * Holwerda et al. (2015) Holwerda, B. W., Bouwens, R., Oesch, P., et al. 2015, ApJ, 808, 6, doi: 10.1088/0004-637X/808/1/6
  * Holwerda et al. (2020) Holwerda, B. W., Bridge, J. S., Steele, R. L., et al. 2020, AJ, 160, 154, doi: 10.3847/1538-3881/aba617
  * Hu et al. (2017) Hu, W., Wang, J., Zheng, Z.-Y., et al. 2017, ApJ, 845, L16, doi: 10.3847/2041-8213/aa8401
  * Inoue et al. (2014) Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805, doi: 10.1093/mnras/stu936
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python
  * Kawamata et al. (2018) Kawamata, R., Ishigaki, M., Shimasaku, K., et al. 2018, ApJ, 855, 4, doi: 10.3847/1538-4357/aaa6cf
  * Kennicutt (1989) Kennicutt, Robert C., J. 1989, ApJ, 344, 685, doi: 10.1086/167834
  * Kennicutt (1998) —. 1998, ApJ, 498, 541, doi: 10.1086/305588
  * Kimble et al. (2008) Kimble, R. A., MacKenty, J. W., O’Connell, R. W., & Townsend, J. A. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7010, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series
  * Kimble et al. (1998) Kimble, R. A., Woodgate, B. E., Bowers, C. W., et al. 1998, ApJ, 492, L83, doi: 10.1086/311102
  * Koekemoer et al. (2003) Koekemoer, A. M., Fruchter, A. S., Hook, R. N., & Hack, W. 2003, in HST Calibration Workshop : Hubble after the Installation of the ACS and the NICMOS Cooling System, 337
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36, doi: 10.1088/0067-0049/197/2/36
  * Kron (1980) Kron, R. G. 1980, ApJS, 43, 305, doi: 10.1086/190669
  * Labbé et al. (2003) Labbé, I., Franx, M., Rudnick, G., et al. 2003, AJ, 125, 1107, doi: 10.1086/346140
  * Laporte et al. (2017) Laporte, N., Nakajima, K., Ellis, R. S., et al. 2017, ApJ, 851, 40, doi: 10.3847/1538-4357/aa96a8
  * Larson et al. (2022) Larson, R. L., Finkelstein, S. L., Hutchison, T. A., et al. 2022, arXiv e-prints, arXiv:2203.08461.  https://arxiv.org/abs/2203.08461
  * Livermore et al. (2017) Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113, doi: 10.3847/1538-4357/835/2/113
  * Livermore et al. (2018) Livermore, R. C., Trenti, M., Bradley, L. D., et al. 2018, ApJ, 861, L17, doi: 10.3847/2041-8213/aacd16
  * Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97, doi: 10.3847/1538-4357/837/1/97
  * Madau (1995) Madau, P. 1995, ApJ, 441, 18, doi: 10.1086/175332
  * Mainali et al. (2018) Mainali, R., Zitrin, A., Stark, D. P., et al. 2018, MNRAS, 479, 1180, doi: 10.1093/mnras/sty1640
  * Mason et al. (2015a) Mason, C. A., Trenti, M., & Treu, T. 2015a, ApJ, 813, 21, doi: 10.1088/0004-637X/813/1/21
  * Mason et al. (2015b) Mason, C. A., Treu, T., Schmidt, K. B., et al. 2015b, ApJ, 805, 79, doi: 10.1088/0004-637X/805/1/79
  * Matthee et al. (2021) Matthee, J., Naidu, R. P., Pezzulli, G., et al. 2021, arXiv e-prints, arXiv:2110.11967.  https://arxiv.org/abs/2110.11967
  * McLure et al. (2010) McLure, R. J., Dunlop, J. S., Cirasuolo, M., et al. 2010, MNRAS, 403, 960, doi: 10.1111/j.1365-2966.2009.16176.x
  * Morishita (2021) Morishita, T. 2021, ApJS, 253, 4, doi: 10.3847/1538-4365/abce67
  * Morishita et al. (2018) Morishita, T., Trenti, M., Stiavelli, M., et al. 2018, ApJ, 867, 150, doi: 10.3847/1538-4357/aae68c
  * Morishita et al. (2020) Morishita, T., Stiavelli, M., Trenti, M., et al. 2020, ApJ, 904, 50, doi: 10.3847/1538-4357/abba83
  * Naidu et al. (2020) Naidu, R. P., Tacchella, S., Mason, C. A., et al. 2020, ApJ, 892, 109, doi: 10.3847/1538-4357/ab7cc9
  * Naidu et al. (2021) Naidu, R. P., Matthee, J., Oesch, P. A., et al. 2021, MNRAS, doi: 10.1093/mnras/stab3601
  * Oesch et al. (2018) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105, doi: 10.3847/1538-4357/aab03f
  * Oesch et al. (2013) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2013, ApJ, 773, 75, doi: 10.1088/0004-637X/773/1/75
  * Oesch et al. (2014) —. 2014, ApJ, 786, 108, doi: 10.1088/0004-637X/786/2/108
  * Oesch et al. (2016) Oesch, P. A., Brammer, G., van Dokkum, P. G., et al. 2016, ApJ, 819, 129, doi: 10.3847/0004-637X/819/2/129
  * Oke & Gunn (1983) Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713, doi: 10.1086/160817
  * Papovich et al. (2016) Papovich, C., Shipley, H. V., Mehrtens, N., et al. 2016, ApJS, 224, 28, doi: 10.3847/0067-0049/224/2/28
  * Patten et al. (2006) Patten, B. M., Stauffer, J. R., Burrows, A., et al. 2006, ApJ, 651, 502, doi: 10.1086/507264
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266, doi: 10.1086/340952
  * Peng et al. (2010) —. 2010, AJ, 139, 2097, doi: 10.1088/0004-6256/139/6/2097
  * Quadri et al. (2007) Quadri, R., Marchesini, D., van Dokkum, P., et al. 2007, AJ, 134, 1103, doi: 10.1086/520330
  * Ravindranath et al. (2020) Ravindranath, S., Monroe, T., Jaskot, A., Ferguson, H. C., & Tumlinson, J. 2020, ApJ, 896, 170, doi: 10.3847/1538-4357/ab91a5
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362, doi: 10.1086/367745
  * Roberts-Borsani et al. (2021) Roberts-Borsani, G., Morishita, T., Treu, T., Leethochawalit, N., & Trenti, M. 2021, arXiv e-prints, arXiv:2106.06544.  https://arxiv.org/abs/2106.06544
  * Rojas-Ruiz et al. (2020) Rojas-Ruiz, S., Finkelstein, S. L., Bagley, M. B., et al. 2020, ApJ, 891, 146, doi: 10.3847/1538-4357/ab7659
  * Ryan et al. (2011) Ryan, R. E., Thorman, P. A., Yan, H., et al. 2011, ApJ, 739, 83, doi: 10.1088/0004-637X/739/2/83
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103, doi: 10.1088/0004-637X/737/2/103
  * Schmidt et al. (2014a) Schmidt, K. B., Treu, T., Trenti, M., et al. 2014a, ApJ, 786, 57, doi: 10.1088/0004-637X/786/1/57
  * Schmidt et al. (2014b) Schmidt, K. B., Treu, T., Brammer, G. B., et al. 2014b, ApJ, 782, L36, doi: 10.1088/2041-8205/782/2/L36
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243, doi: 10.1086/146614
  * Schmidt (1963) —. 1963, ApJ, 137, 758, doi: 10.1086/147553
  * Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65, doi: 10.1086/373922
  * Somerville et al. (2008) Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481, doi: 10.1111/j.1365-2966.2008.13805.x
  * Somerville et al. (2015) Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337, doi: 10.1093/mnras/stv1877
  * Stark (2016) Stark, D. P. 2016, ARA&A, 54, 761, doi: 10.1146/annurev-astro-081915-023417
  * Stefanon et al. (2021) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, arXiv e-prints, arXiv:2103.16571.  https://arxiv.org/abs/2103.16571
  * Stefanon et al. (2019) Stefanon, M., Labbé, I., Bouwens, R. J., et al. 2019, ApJ, 883, 99, doi: 10.3847/1538-4357/ab3792
  * Steidel et al. (2016) Steidel, C. C., Strom, A. L., Pettini, M., et al. 2016, ApJ, 826, 159, doi: 10.3847/0004-637X/826/2/159
  * Tacchella et al. (2022) Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170, doi: 10.3847/1538-4357/ac4cad
  * Tang et al. (2018) Tang, M., Stark, D., Chevallard, J., & Charlot, S. 2018, ArXiv e-prints.  https://arxiv.org/abs/1809.09637
  * Tilvi et al. (2016) Tilvi, V., Pirzkal, N., Malhotra, S., et al. 2016, ApJ, 827, L14, doi: 10.3847/2041-8205/827/1/L14
  * Trenti & Stiavelli (2008) Trenti, M., & Stiavelli, M. 2008, ApJ, 676, 767
  * Trenti et al. (2011) Trenti, M., Bradley, L. D., Stiavelli, M., et al. 2011, ApJ, 727, L39, doi: 10.1088/2041-8205/727/2/L39
  * Trenti et al. (2012) —. 2012, ApJ, 746, 55, doi: 10.1088/0004-637X/746/1/55
  * Vogelsberger et al. (2020) Vogelsberger, M., Nelson, D., Pillepich, A., et al. 2020, MNRAS, 492, 5167, doi: 10.1093/mnras/staa137
  * Whitaker et al. (2011) Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86, doi: 10.1088/0004-637X/735/2/86
  * Wilkins et al. (2017) Wilkins, S. M., Feng, Y., Di Matteo, T., et al. 2017, MNRAS, 469, 2517
  * Yan et al. (2011) Yan, H., Yan, L., Zamojski, M. A., et al. 2011, ApJ, 728, L22
  * Yung et al. (2019a) Yung, L. Y. A., Somerville, R. S., Finkelstein, S. L., Popping, G., & Davé, R. 2019a, MNRAS, 483, 2983, doi: 10.1093/mnras/sty3241
  * Yung et al. (2019b) Yung, L. Y. A., Somerville, R. S., Popping, G., et al. 2019b, arXiv e-prints, arXiv:1901.05964.  https://arxiv.org/abs/1901.05964
  * Zitrin et al. (2015) Zitrin, A., Labbé, I., Belli, S., et al. 2015, ApJ, 810, L12, doi: 10.1088/2041-8205/810/1/L12

##  Appendix A Candidate Photometry and Luminosity Function Measurements

In Table 6, we provide the photometry of each candidate in all available
filters in units of nJy. In Table 7, we report the median and 68% range of the
number density posterior distribution plotted in Figure 17.

Table 6: Sample Photometry

=5mm  Par ID $V_{606}$ $Y_{098}$ $J_{125}$ $H_{160}$ 3.6µm 4.5µm (nJy) (nJy)
(nJy) (nJy) (nJy) (nJy) Par1033+5051 116 $-0.35\pm 4.71$ $0.11\pm 9.47$
$153.80\pm 11.21$ $145.76\pm 15.83$ $420.31\pm 188.84$ $1386.37\pm 744.31$
Par0456-2203 473 $-26.35\pm 21.49$ $-10.55\pm 16.58$ $176.79\pm 31.70$
$563.72\pm 50.93$ $2964.83\pm 668.70$ $2987\pm 1079.72$ Par0926+4000 369
$23.34\pm 27.31$ $-6.60\pm 11.81$ $51.28\pm 36.70$ $487.13\pm 35.67$
$162.72\pm 356.68$ $\cdots$ Par0756+3043 92 $-1.17\pm 12.04$ $1.22\pm 15.25$
$-9.08\pm 11.69$ $140.53\pm 25.49$ $276.56\pm 455.82$ $389.83\pm 1186.16$ Par
ID $I_{600}$ $Y_{098}$ $J_{125}$ $H_{160}$ 3.6µm 4.5µm Par0440-5244 497
$-9.97\pm 18.86$ $-6.55\pm 9.49$ $193.34\pm 18.95$ $235.66\pm 29.11$
$333.83\pm 1767.52$ $738.38\pm 623.41$ Par1301+0000 37 $-24.32\pm 43.33$
$-3.41\pm 16.46$ $-5.51\pm 24.00$ $246.95\pm 33.66$ $-340.46\pm 437.01$
$1542.47\pm 1349.34$ Par0440-5244 593 $-2.27\pm 7.16$ $7.23\pm 5.21$ $-0.97\pm
4.72$ $94.47\pm 10.81$ $153.45\pm 290.15$ $-318.25\pm 1150.19$ Par ID
$I_{600}$ $Y_{105}$ $J_{125}$ $H_{160}$ 3.6µm 4.5µm Par0713+7405 95 $38.80\pm
28.05$ $-3.33\pm 14.72$ $107.11\pm 21.64$ $183.20\pm 23.86$ $251.30\pm 217.45$
$277.66\pm 661.93$ Par0259+0032 194 $-3.20\pm 20.21$ $-9.25\pm 12.88$
$148.17\pm 23.77$aaThe $J_{125}$ photometry quoted here for Par0259+0032_194
is that measured after masking out the bad pixels as described in Section
4.1.5. $440.55\pm 40.74$ $\cdots$ $\cdots$ Par0843+4114 120 $-20.46\pm 13.94$
$21.41\pm 16.82$ $-1.83\pm 9.23$ $101.05\pm 19.05$ $-190.14\pm 521.57$
$-138.04\pm 1263.79$ Par0956-0450 684 $-6.77\pm 17.48$ $-8.50\pm 11.00$
$10.04\pm 12.52$ $404.85\pm 18.50$ $16.36\pm 481.12$ $335.27\pm 707.27$
Par0926+4536 155 $-4.68\pm 12.92$ $15.20\pm 15.14$ $-16.90\pm 8.77$ $102.65\pm
18.64$ $144.67\pm 159.01$ $240.19\pm 746.29$ Par ID $V_{606}$ $I_{814}$
$J_{110}$ $H_{160}$ 3.6µm 4.5µm Par335 251 $\cdots$ $-1.64\pm 14.63$ $14.54\pm
17.31$ $460.15\pm 34.90$ $77.47\pm 719.45$ $\cdots$

Note. — The candidates are ordered by filter set, eather than in redshift-
order as in Table 5.

Table 7: Measured Number Densities for $z\sim 8.5-11$

=5mm  UV Magnitude Number Density 68% Range ($10^{-6}$ Mpc-3) ($10^{-6}$
Mpc-3) -23.15 $<0.91$ $\cdots$ -23.05 1.02 $0.68-1.49$ -22.95 1.18 $0.77-1.72$
-22.85 1.33 $0.88-1.95$ -22.75 1.44 $0.96-2.10$ -22.65 1.48 $0.97-2.16$ -22.55
1.54 $1.00-2.27$ -22.45 1.66 $1.08-2.49$ -22.35 1.92 $1.25-2.85$ -22.25 2.21
$1.41-3.27$ -22.15 2.53 $1.63-3.76$ -22.05 2.93 $1.91-4.37$ -21.95 3.48
$2.26-5.16$ -21.85 4.22 $2.74-6.32$ -21.75 5.17 $3.33-7.71$ -21.65 6.35
$4.09-9.46$ -21.55 7.86 $5.07-11.70$ -21.45 9.98 $6.48-14.83$ -21.35 12.98
$8.37-19.12$ -21.25 16.87 $11.00-24.73$ -21.15 20.40 $13.16-30.36$ -21.05
23.44 $15.17-34.63$ -20.95 25.95 $16.75-38.65$ -20.85 28.98 $18.64-43.52$
-20.75 37.59 $23.70-56.79$ -20.65 50.39 $31.65-76.55$

Note. — The number densities are taken as the median of the MCMC posterior at
each magnitude step using the “pseudo-binning” technique of Finkelstein et al.
(2022). The 68% range of the posterior distribution listed in the third column
is plotted as the dark blue shaded region in Figure 17.

##  Appendix B Rejected candidates

Figure 18: The $H_{160}$ postage stamps of the 97 spurious sources rejected in
the non-subjective categories during our source vetting (Sections 3.4 and
3.3). Group (a) shows $5\arcsec\times 5\arcsec$ stamps of the 48 sources that
were rejected as satellite trail remnants or diffraction spikes through visual
inspection. Group (b) shows $5\arcsec\times 5\arcsec$ stamps of the 12 sources
that were rejected along image edges or near image edges in regions of
increased noise. Group (c) shows stamps of the 25 sources rejected as either
bad pixels or sources with photometry contaminated by bad pixels through
visual inspection of the RMS maps at the source positions. We adpot a
threshold of $RMS>0.1$, equivalent to weights of $<100$. The bad pixels are
identified by red contours. Each stamp is 2″ on a side, zooming in on the
source position to highlight the source shapes. Group (d) shows
$5\arcsec\times 5\arcsec$ stamps of the 12 sources that were rejected as image
persistence. Ten were identified based on analysis of counts in FLTs observed
in the 24 hours preceding the observation of the FLTs that contributed to the
mosaic. The remaining two (surrounded by red boxes) were identified in the
WISP survey based on the source position relative to a bright source. Each
stamp is displayed on a zscale interval (calculated using the pixels in the 5″
stamp) with no image smoothing or pixel interpolation. The $r=0\farcs 75$
circles show the position of the rejected candidate as a reference.

In this section, we show $H_{160}$ postage stamps for all 166 rejected
sources. These candidates were selected via the criteria discussed in Section
3.2 and then later rejected through persistence checks (Section 3.3), or
visual inspection (Section 3.4). All stamps are displayed on a linear stretch
with a zscale interval that highlights pixel values near the image median.
These stamp normalizations are calculated using the pixels in a
$5\arcsec\times 5\arcsec$ box centered on the source position. We note that on
the reverse colormap adopted here, very negative values are displayed as white
pixels. We do not apply any pixel interpolation or smoothing to these stamps,
and use $r=0\farcs 75$ circles to indicate the source positions.

Of the 193 candidates that pass our selection criteria, the majority of them
are spurious detections that we classify visually in four categories. First,
we reject 48 sources as satellite remnants or diffraction spikes (Figure 18a).
As can be seen in Figure 18a, many of the sources in this category are
detected as part of the same set of satellite trails. Next, we reject 12
sources that are detected on or very close to an image edge (Figure 18b). The
sources in this category suffer from incomplete source extraction and
unreliable photometry. Through inspection of the RMS maps at the position of
each candidate, we identify 25 sources that are bad pixels in the RMS map
(RMS$>$0.1, equivalent to a weight less than 100) or are contaminated by bad
pixels (i.e., bad pixels are included in the isophotal area of the source,
therefore impacting the photometry). We show the central $2\arcsec\times
2\arcsec$ around each of these candidates in Figure 18c to zoom in on the
morphologies of these spurious or contaminated sources, and indicate the
locations of the bad pixels with red contours. The final category of visually
rejected sources are those that have been identified as hot pixels, detector
artifacts, or sources with strange morphologies. This category is the only
subjective part of our visual inspection, and we attempt to account for these
rejections in our calculate of the effective volume (see Section 6.3). We show
the 81 sources rejected in this category in Figure 19, zooming in on the
celtral $2\arcsec\times 2\arcsec$ to highlight the morphologies.

Finally, we reject 12 candidates as contamination by image persistence. In
these cases, a previous observation of a bright source left a residual charge
on the detector that fades with time. If this leftover flux is detected in the
$H_{160}$ and is absent from the other filters, the persistence will be
selected as a high-redshift source. We have identified 10 cases of persistence
based on the counts registered on the detector at the position of the
candidates prior to the candidate’s observation. An additional 2 WISP sources
are identified as persistence due to their location relative to a bright
source and the survey strategy of the WISP observations (see Section 3.3). We
show all 12 cases of persistence in Figure 18d.

Figure 19: The $H_{160}$ postage stamps ($2\arcsec\times 2\arcsec$) of the 81
sources that were rejected as having unreliable source shapes or morphologies.
We note that this rejection is the only subjective part of the visual
inspection and rejection process, and we quantify the biases introduced
through these rejections using Figure 15. The stamps are displayed in the same
manner as those in Figure 18. The two candidates from B16 that we rejected
during visual inspection are identified by red boxes: borg_0240-1857_25 in the
top row (third from the right) and borg_0456-2203_1091 in the second row
(fifth from the left).