# SKYSURF-3: Testing Crowded Object Catalogs in the Hubble eXtreme Deep Field
Mosaics to Study Sample Incompleteness from an Extragalactic Background Light
Perspective

Darby M. Kramer School of Earth and Space Exploration, Arizona State
University, Tempe, AZ 85287 USA  Timothy Carleton School of Earth and Space
Exploration, Arizona State University, Tempe, AZ 85287 USA  Seth. H. Cohen
School of Earth and Space Exploration, Arizona State University, Tempe, AZ
85287 USA  Rolf A. Jansen School of Earth and Space Exploration, Arizona State
University, Tempe, AZ 85287 USA  Rogier A. Windhorst School of Earth and Space
Exploration, Arizona State University, Tempe, AZ 85287 USA  Norman A. Grogin
Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218,
USA  Anton M. Koekemoer Space Telescope Science Institute, 3700 San Martin
Dr., Baltimore, MD 21218, USA  John W. Mackenty Space Telescope Science
Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA  Norbert Pirzkal
Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218,
USA

###### Abstract

Extragalactic Background Light (EBL) studies have revealed a significant
discrepancy between direct measurements — via instruments measuring “bare” sky
from which Zodiacal and Galactic light models are subtracted — and
measurements of the Integrated Galaxy Light (IGL). This discrepancy could lie
in either method, whether it be an incomplete Zodiacal model or missed faint
galaxies in the IGL calculations. It has been proposed that the discrepancy is
due to deep galaxy surveys, such as those with the _Hubble Space Telescope_
(_HST_), missing up to half of the faint galaxies with 24 $\sim$$<$ $m_{\rm
AB}$ $\sim$$<$ 29 mag. We address this possibility by simulating higher number
densities of galaxies, and so assess incompleteness due to object overlap,
with three replications of the Hubble UltraDeep Field (HUDF). SourceExtractor
is used to compare the recovered counts and photometry to the original HUDF,
allowing us to assess how many galaxies may have been missed due to confusion,
i.e., due to blending with neighboring faint galaxies. This exercise reveals
that, while up to 50% of faint galaxies with 28 $\sim$$<$ $m_{\rm AB}$
$\sim$$<$ 29 mag were missed or blended with neighboring objects in certain
filters, not enough were missed to account for the EBL discrepancy alone in
any of the replications.

Cosmology: Extragalactic Background Light — Galaxies: Galaxy Counts — Solar
System: Zodiacal Light — Instrument: Hubble Space Telescope

††software: SourceExtractor (Bertin & Arnouts, 1996)

##  1 Introduction

Extragalactic Background Light (EBL) is the dominant background in the
Universe after the cosmic microwave background (CMB), and is generally
considered to be comprised of the short wavelength (0.1 µm) ultraviolet to
long wavelength (1000 µm) far-infrared emission (Driver et al., 2016; Mattila
& Väisänen, 2019), though this depends on convention, and can also be
considered to include the entire electromagnetic spectrum. In terms of energy
density, two main components make up the EBL: the cosmic optical background
(COB) and the cosmic infrared background (CIB). Because a range of
astrophysical processes emit significant amounts of photons at ultraviolet
(UV), optical and infrared (IR) wavelengths, the EBL includes radiation from
several types of objects in the Universe such as stars, AGN, and dust (Madau &
Pozzetti, 2000). There are two main ways that EBL measurements are obtained:
(A) direct measurements of the sky, from which Zodiacal and diffuse Galactic
light is subtracted, and (B) integrated galactic light (IGL) from galaxy
surveys. Both methods are demonstrated in, e.g., Bernstein (2007). The EBL is
expected to be highly dominated by galaxy light, and current technology allows
for deep images at UV, optical, and IR wavelengths. Thus, integrating the
galaxy light in deep sky catalogs should yield an accurate representation of
the EBL, unless a significant diffuse EBL component is not present in the
discrete galaxy catalogs.

The reason this is still an active area of study is that, although the IGL at
IR wavelengths is in fairly good agreement with direct CIB measurements, the
UV and optical IGL still do not agree with the directly-estimated diffuse COB:
current direct measurements of the EBL are significantly higher than the IGL
calculations and COB models (Keenan et al., 2010; Driver et al., 2016; Hill et
al., 2018). Uncertainties in the direct EBL measurements come from foregrounds
such as Zodiacal light and dust from the Milky Way (Madau & Pozzetti, 2000;
Matsumoto et al., 2015), while uncertainties in IGL calculations come from
instrumental limitations and the number density of galaxies on the sky.

Driver (2021) addresses the issue with high energy gamma ray studies, such as
HESS and MAGIC, which get information from blazar photons that interact with
the EBL on their path toward Earth. The blazar spectra are expected to have
absorption features at the wavelengths where the gamma ray and EBL photons
interact, constraining the permitted levels of the EBL intensity with a
$\sim$50% error range (Aharonian et al., 2006; Acciari et al., 2019). These
predictions match the IGL calculations and EBL models well (Abeysekara et al.,
2019), providing independent constraints that suggest the higher direct EBL
measurements may still include components that need to be subtracted. Deep
space missions, such as Pioneer 10/11 and New Horizons, also improve direct
EBL measurements because they have much less Zodiacal light to contend with
than telescopes closer to the Sun (Matsuoka et al., 2011; Lauer et al., 2021,
2022). These missions have provided some direct EBL data points which are
lower than those from Low Earth Orbit or ground-based experiments (see, e.g.,
Bernstein, 2007) Even with these improved constraints, the cause of the EBL
discrepancy is still unidentified (Driver, 2021). Among the possibilities are:
(A) a Zodiacal component that is missing from or underestimated in current
models, causing an overestimation of the EBL in direct measurements; (B)
missed faint diffuse galaxies or extended galaxy outskirts in galaxy surveys,
causing an underestimation of the IGL; or (C) both.

Lauer et al. (2021) propose that their extra, unexplained diffuse flux
component of the COB might result if galaxy surveys were missing up to half of
all the faint galaxies that exist. If these faint galaxies are being missed,
they are not counted in the IGL calculations and could cause the true EBL to
be underestimated. Their paper explains that the faint-end slope could drop
off sharply at $m_{\rm AB}$ $>$ 24 mag, i.e., that surveys miss significant
numbers of faint galaxies with $m_{\rm AB}$ $\sim$$<$ 30 mag. Such galaxies
are very faint, however, and contribute much less than 25% of the IGL
according to Windhorst et al. (2022). That paper and Carleton et al. (2022)
suggest that one would need a factor of 4–8 times more faint galaxies with 24
$\sim$$<$ $m_{\rm AB}$ $\sim$$<$ 30 mag to make up for the missing diffuse EBL
flux (Lauer et al., 2021, 2022; Bernstein, 2007; Hill et al., 2018).

Deep extragalactic surveys with _HST_ , starting with the Hubble Deep Field
(Williams et al., 1996), and subsequently the Hubble UltraDeep Field/eXtreme
Deep Field (HUDF/XDF) (Beckwith et al., 2006; Oesch et al., 2010; Ellis et
al., 2013; Koekemoer et al., 2013; Illingworth et al., 2013; Rafelski et al.,
2015) have detected and resolved distant galaxies, with some candidates above
z $\sim$ 10 and to apparent magnitudes $m_{\rm AB}$ $\sim$ 30 mag,
constraining the observed galaxy luminosity function (Oesch et al., 2009;
Finkelstein et al., 2012; Windhorst et al., 2011, 2021). Even though new
galaxy surveys continue to reach deeper and wider on the sky, their galaxy
counts are still complicated by instrumental effects and object confusion.
Instrumental confusion is the inability of the optical system of telescope,
camera, and detector to resolve and separate faint objects because of limits
that depend on the image point-source sensitivity. On the contrary, natural
confusion is the inability to resolve faint, extended objects due to their
surface brightness and statistical overlap. Any faint object that is even
slightly resolved can overlap other objects. This causes incompleteness at dim
magnitudes.

Although instrumental confusion can be mitigated with larger apertures,
natural confusion is a fundamental problem that will persist as galaxy surveys
reach deeper. Confusion is thus a major source of noise in the deep images
being acquired today (Roseboom et al., 2010). This work focuses on generating
and studying the effects of natural confusion alone, as it has stronger limits
on present-day surveys. Deep surveys with _HST_ run into the natural confusion
limit of 1 object per 25 – 50 independent beams, where the beam is the average
object diameter at a given flux (Serjeant et al., 1997; Silva et al., 2005;
Windhorst et al., 2008, their Figure 3).

Figure 1: Examples of square image regions excised for analysis from the XDF
mosaics before the objects brighter than $m_{\rm AB}$ $\sim$ 24 mag were
masked out. (_a_) Exposure times in the WFC3/IR images, such as the F160W
image shown here, vary significantly throughout a single mosaic. We therefore
excised a large region with relatively uniform depth in the deepest part of
the WFC3/IR mosaics — a 1401 $\times$ 1401 pixel region centered on pixel
(2734,3072) in the XDF mosaics. (_b_) We similarly excised the largest square
image regions with relatively uniform depth from the ACS/WFC mosaics, such as
the F606W image shown here. Each $\sim$ 2870 $\times$ 2870 pixel region is
centered on pixel (2625,2625) in the XDF mosaics.

This is the third paper associated with the SKYSURF project, an Archival _HST_
program (AR 15810; PI: R. Windhorst) that aims to resolve the direct EBL vs.
IGL conflict. SKYSURF is addressing the conflict by re-processing _HST_ data
in order to measure the sky values, instead of reducing them (Windhorst et
al., 2022; Carleton et al., 2022).

Windhorst et al. (2022) showed that bright objects with $m_{\rm AB}$$<17$ mag
contribute $\sim 25\%$ of the IGL, while objects 17 $\sim$$<$ $m_{\rm AB}$
$\sim$$<$ 22 mag contribute about 50% and objects fainter than $m_{\rm
AB}$$=22$ mag contribute 25% of the IGL. Moreover, objects fainter than
$m_{\rm AB}$$=24$ contribute less than 10% of the IGL. Therefore, in order to
explain a discrepancy as large as a factor of 3–5 between the direct estimates
of the EBL and estimates of the EBL from IGL measurements, one would require
as many as 8 times the number of faint galaxies ($m_{\rm AB}$$>24$ mag) than
are currently detected.

In this work, we study the XDF images that were released by Illingworth et al.
(2013) to explore how object confusion could affect estimates of the
extrapolated galaxy light. We investigate specifically how many times we can
replicate the XDF data onto itself, after successive 90∘ rotations, before
object overlap starts to significantly limit the ability of the deblending
algorithm in SourceExtractor to construct complete object catalogs. In § 2 we
discuss the _HST_ data used in this work and how it was pre-processed. § 3
discusses the use of SourceExtractor on these images. In § 4 we show and
analyze the SourceExtractor source photometry extracted from the deep field
images, and in § 5 the conclusions and implications of these results are laid
out.

##  2 Data

It was desired to perform this exercise with data rather than simulations
(such as those in Windhorst et al. (2022) § 4.2) so that galaxy sizes and
morphologies could be as realistic as possible. We therefore use the XDF
dataset that was released by Illingworth et al. (2013). These images contain
data obtained with the Advanced Camera for Surveys Wide Field Channel
(ACS/WFC, Ryon, 2022) between 2002 and 2013, including new and reprocessed ACS
data obtained after the original HUDF to improve its depth, as well as data
obtained with the Wide Field Camera 3 InfraRed Channel (WFC3/IR, Dressel,
2022) between 2009 and 2013. The XDF ACS mosaics are deeper by $\sim 0.12-0.2$
magnitudes compared to the original HUDF (Beckwith et al., 2006), while the
XDF WFC3/IR mosaics reach similar depths to the mosaics released from the
HUDF12 campaign (Ellis et al., 2013; Koekemoer et al., 2013) but were reduced
completely independently, with both of these releases containing the full set
of WFC3/IR data obtained on this field at that time. For consistency with the
XDF ACS mosaics, we also make use of the XDF WFC3/IR mosaics in this work. We
start with the fully-reduced public XDF mosaics in the F435W, F606W, and F775W
broad-band filters, but more filters were used in the analysis (F105W, F125W,
F160W). The F140W, F814W, and F850LP filters were omitted from this study
because they are shallower than the others used here. The _HST_ data used in
this work can be found in MAST: http://dx.doi.org/10.17909/T97P46 (catalog
10.17909/T97P46).

To avoid unnecessary complications from the severely declining stacked
exposure depth around the edges of the field, we excised image regions from
all filters’ mosaics that were exactly square and had near-uniform exposure
depth. The WFC3/IR mosaics were masked to only include a region with
relatively uniform depth in the deepest part of the F105W, F125W, and F160W
images — a 1401$\times$1401 pixel region centered on pixel (2734, 3072). Fig.
1 _a_ shows the section that was used compared to the full XDF mosaic. The
same square image regions were excised from the associated weight images.

The ACS/WFC F435W, F606W, and F775W mosaics were masked to excise squares
defined to have diagonals measuring 4059 pixels, making their sides
approximately 2870 pixels in length (Fig. 1 _b_). Each excised $\sim
2870\times 2870$ pixel region is centered on pixel (2625, 2625) in the
original XDF mosaics.

Figure 2: A centered, zoomed comparison of the four images used for the F606W
analysis. Each detected object is marked with a green ellipse, the size of
which corresponds to $2.5\times$ the profile rms along major and minor axes to
account for the Kron factor (Kron, 1980). (_a_) Original XDF image with
objects brighter than $m_{\rm AB}$ $\sim$ 24 mag masked out. (_b_) The masked
original image + a masked and _once_ rotated original image. (_c_) The image
of panel (b) + a masked and _twice_ rotated original image. (_d_) The image of
panel (c) + a masked and _three times_ rotated original image. This process
increases the surface density of objects brighter than the detection limit by
integer factors in each iteration, while leaving the sky unchanged.

##  3 Methods

SourceExtractor (Bertin & Arnouts, 1996) has been used by astronomers for
decades to identify objects in images and extract relevant object properties,
including object photometry. We used SourceExtractor along with Python to
process and analyze each image following established methods for taking
associated weight images into account, and we will use MAG_AUTO object
magnitudes throughout. In this work, we define “sky” or “sky values” to mean
pixels that were not assigned to part of an object by SourceExtractor.

We wish to investigate only galaxies fainter than $m_{\rm AB}$ = 24 mag
because of the possibility that significant numbers of objects in this range
are missed due to confusion. For this reason, the objects brighter than
$m_{\rm AB}$ = 24 mag were masked from the original XDF data before rotation
and addition. To do this, Gaussian distributions were constructed with the sky
and root-mean-square (rms) values returned by SourceExtractor for each filter.
Because, for example, the 174.4 ks F606W image consists of 286 separate
exposures, its noise and weight map are very uniform. Therefore, we drew
random sky + noise values from our constructed distribution and placed these
in the pixels where the bright objects had been identified by a segmentation
map and removed. The segmentation map was created with detection and analysis
thresholds of $0.8\times$ the sky rms value (which is $\sim$53% lower than the
$1.5\sigma$ detection values over 4 connected pixels used in the analysis of
the final images here) in order to detect and also mask the extended light
surrounding bright galaxies. This method eliminated the strong object overlap
from large, bright galaxies, and the result is shown in Figure 2 for the
masked F606W image and replications thereof. In all, between 3 and 10% of the
pixels in these mosaics were masked.

To prepare for rotation and addition, the original images were run through
SourceExtractor with the option to return FITS files with only the objects in
them. These object-only FITS images were rotated and co-added to the original
mosaics. This minimizes the negative effects of co-adding sky noise.

When the ACS/WFC and WFC3/IR images were drizzled from their native pixels to
the 60 milliarcsecond square pixel images (used in this work), the resultant
pixels had necessarily correlated noise (Casertano et al., 2000). Moreover,
image rotation by any angle other than $n$ $\times$ 90∘ and stacking would
introduce additional noise due to interpolation between unaligned pixels,
whether square or not. We therefore only used image transpositions and
inversions to implement the rotations over $n$ $\times$ 90∘. Each “rotated”
image was co-added to the original. This multiplied the number of potentially
detectable objects by 2$\times$, 3$\times$, or 4$\times$, respectively.
Because the weight maps are inverse variance weights, the sum of the inverses
of the original weight maps was inverted to get the replication weight maps.
That is, for the 3$\times$ replication, the weight image was constructed from
the original XDF weight map and the rotated weight maps:

$\text{W3x}=\frac{1}{1/\text{W}_{0}+1/\text{W}_{90}+1/\text{W}_{180}}\ ,$ (1)

where $\text{W}_{90}$ is the weight image corresponding to a 90∘ rotation of
the original XDF weight image. SourceExtractor then processed all the images
with their appropriate zero points and weight images found on the XDF Data
Release page; detection and analysis thresholds set to $1.5\sigma$; the
minimum number of connected pixels above these thresholds to be a detection
set to 4; and minimum contrast for deblending set to 0.02.

SourceExtractor counted the objects in the masked original XDF data, as well
as in the replicated data, and returned object catalogs for each. The
magnitude and half-light radius columns of the catalog files were used to
construct the following figures. Elliptical apertures were generated using the
position, size, and orientation columns for each detected object, and imported
into DS9 (Joye & Mandel, 2003) for visual inspection of the objects detected
in the respective images. Figure 2 shows examples of these for the F606W
filter for all four realizations.

Error calculations in this analysis were performed as follows. For the
histograms (left panels in Fig. 3), the error bars were calculated simply by
$\sigma_{\rm bin}=\sqrt{n}$, where $n$ is the number of objects in a
particular bin of the original image. For the ratios of the histograms (the
right panels in Fig. 3), the uncertainties are based on their corresponding
histograms and derived for a Gaussian error distribution:

$\sigma_{\rm
bin}=\frac{m}{pn}*\sqrt{\left(\frac{\sqrt{m}}{m}\right)^{2}+\left(\frac{\sqrt{pn}}{pn}\right)^{2}}$
(2)

where $p$ is the integer by which the original image bins are multiplied (the
effective multiple in the transposed and added images), and $m$ is the number
of objects in a particular bin of the rotated and co-added images. We note
that the rms noise increased by 4–7% from the initial XDF mosaics to their
$4\times$ replications. This can be explained by the wings of some objects
going undetected by SourceExtractor, and therefore being included in sky
statistics.

We also tested whether using translation rather than rotation affects the
results of an analysis like the present one. It turns out that, although the
rotation causes slight statistical distortion at the image center, the
difference in the results between the two methods is within the counting
errors of $\sqrt{n}$ for each bin, and generally well-within 10%. As we will
see below, this difference between the two replication approaches is much less
than the incompleteness levels caused by object overlap that we seek to
quantify below. Therefore, we chose to maintain the rotation and addition
method, because it dealt well with edge effects (i.e., without cutting
additional galaxies in half).

##  4 Results

We now wish to examine if the transposed and co-added images indeed contain
2$\times$, 3$\times$, or 4$\times$ as many objects as the original,
respectively.

Figure 3: _Left Column:_ Histograms of object counts per $\frac{1}{2}$ mag bin
per deg2 for the F435W, F606W, and F775W original and replicated images up to
their completeness limits. AX Rep refers to the integer A replications of the
XDF data, and AxXDF refers to that integer times the original XDF bin counts.
Points are located at the centers of their bins. Spline fits show the general
trends, and counts are plotted and listed up to their faintest complete
magnitudes for the individual filters (e.g., up to $m_{\rm AB}$ = 29 mag for
F606W). _Right Column:_ The ratios per $\frac{1}{2}$ mag bin of the replicated
image bin counts over two, three, and four times the original image bin
counts. These curves demonstrate the completeness of each replication. The
faint horizontal blue line marks the expected ratio of 1, as a comparison to
the actual values.

Figure 2 allows a visual comparison of these four images for the F606W filter.
The images include (as green ellipses) the apertures associated with each of
the objects detected and measured by SourceExtractor in each of the four
images. Figure 3 shows the results from the above methods as applied to XDF
F435W, F606W, and F775W images in the HUDF. The left column of Figure 3 shows
the differential object counts as a function of $m_{\rm AB}$ for each of the
four realizations, as indicated by the four different colors. These are the
raw data before correcting for the fact that each of the realizations has
2$\times$, 3$\times$, or 4$\times$ as many objects as the original image
(shown in black). Some objects sufficiently overlap neighbors such that
SourceExtractor deemed them to be a single object, causing the count slopes to
slightly change between each iteration. While the (black) input catalog had no
objects brighter than $m_{\rm AB}$ = 24 mag, object overlap even resulted in
some objects being slightly brighter than this limit (these points are present
but excluded in the figures here). The panels in the right column of Figure 3
show the ratios of the counts in the left panels and the original counts (in
black), but multiplied by the proper integer. Therefore, these “completeness
functions” in the right panels show the decreasing trend in completeness as
each replication becomes more and more dense.

Figure 4 shows a magnitude-radius plot for all four F606W samples, which
allows us to compare these samples according to their angular size and
magnitude. The nearly horizontal limit that ends at $m_{\rm AB}$ $\simeq$ 30.5
mag is the point-source sensitivity limit. The slanted cyan line that extends
into the upper right of Figure 4 is the surface brightness limit that
indicates how well larger but dim objects can be resolved in the F606W image
with its 174.4 ks total exposure time. The slanted pink band indicates the
approximate natural confusion limit of 1 object / 25 beams to 1 object / 50
beams, to the right of which objects cannot be easily detected because they
start to significantly overlap in the SourceExtractor segmentation maps (e.g.,
Windhorst et al., 2008). Note that in the concept of natural confusion of
resolved objects, “beam” refers to the average object area at that AB-
magnitude level, or approximately $\pi r_{e}^{2}$, where $r_{e}$ is the median
object half-light radius at that flux level in Figure 4.

Figure 4: Magnitude vs. Half-Light Radius plot for the F606W original and
replicated images. The black points are the original objects, and the other
colors are the 2$\times$, 3$\times$, and 4$\times$ XDF replications’ objects.
Each successive color has more transparent markers, and there are more blue
marks than green, more green than red, and more red than black. This figure
was modeled after Figure 3 in Windhorst et al. (2008) so features could be
marked and compared. The pink band represents the natural confusion limit with
a range between 1/25 and 1/50 objects per beam (see e.g. Serjeant et al.,
1997; Silva et al., 2005). The cyan dashed curve represents the 592-hour
surface brightness limit of XDF data. Its shape has been changed from
Windhorst et al. (2008) as the curve is only intended to draw the eye to the
trend and the cutoff of the objects, demonstrating the exposure time and
resolution limits to observing these objects. The dark grey dashed line
represents the limit below which point sources with fainter magnitudes cannot
be resolved, as they blend in with the sky distribution. All curves have been
shifted in magnitude to account for the $m_{Vega}$ to $m_{\rm AB}$ conversion.
This figure demonstrates that natural confusion (pink band) is an important
completeness constraint for bright, large objects that make up most of the IGL
($m_{\rm AB}$ $\sim$$<$ 28.5 mag). It also demonstrates that the three main
limits on the original HUDF sample remain the the same for the replicated
samples, validating our methods of replication, namely the point source
detection limit (horizontal grey dashed line), the SB-limit (slanted cyan
dashed line) and the natural confusion limit due to statistical object overlap
(pink band).

=-0.55in  Sample Filter F435W⋆ F606W⋆ F775W⋆ F105W† F125W† F160W† HUDF09 0.237
0.228 0.222 0.209 0.226 0.209 XDF 0.241 0.248 0.238 0.249 0.249 0.232
2$\times$XDF 0.223 0.230 0.218 0.211 0.219 0.206 3$\times$XDF 0.213 0.213
0.203 0.178 0.191 0.163 4$\times$XDF 0.206 0.190 0.191 0.149 0.173 0.129 ⋆
indicates $m_{\rm AB}$ 25-29 slope † indicates $m_{\rm AB}$ 25-28 slope

Table 1: Slopes (in dex per 0.5 magnitude) of the samples in this work as
compared to those measured by Windhorst et al. (2011) in the HUDF09 field
(Bunker et al., 2010). XDF count slopes in this work are generally larger than
those in Windhorst et al. (2011) Figure 12, which is likely due to a higher
percentage of completeness in the XDF data. The replicated slopes in this work
are also expected to be flatter with each successive replication of the XDF
data, as is indeed observed.

Comparing slopes in the left column of Figure 3 to Figure 11 of Windhorst et
al. (2011), we find that those in the present work are generally flatter (see
Table 1). Because only fainter magnitudes were used here, this is expected.
First, the slopes are expected to flatten due to incompleteness at the fainter
magnitudes. Second, the same flattening trend is also seen between iterations
of replicating the images due to the increasing severity of the simulated
natural confusion resulting in increasing incompleteness. For the
corresponding results for the WFC3/IR F105W, F125W, and F160W filters, we
refer the reader to the Appendix and Table 1. These results are consistent
with the findings already presented, and all six filters will be discussed in
§ 5.

It is worth noting that the exercises in this paper were initially performed
on un-masked XDF data. Masking brighter objects was later desired to narrow
the focus of this work toward the faint objects that are of importance in this
discussion. Masking the brighter objects before rotation and addition and
comparing to un-masked replications showed that the un-masked analysis had
$\sim$5-30% lower total relative number counts (right column of Figures 3 and
5), meaning fewer objects were recovered in the un-masked images than in the
masked ones. When these object counts were weighted by their luminosities and
summed to get faint-end IGL calculations, the masked IGL values for each
filter for objects with 24 $\sim$$<$ $m_{\rm AB}$ $\sim$$<$ 29 mag were on
average 10% lower than those from the un-masked images. This is due to the
generally-larger number of objects in the un-masked image. Therefore, our
completeness curves and summed values in each _HST_ filter, in the right
panels of Figures 3 and 5, are upper limits.

##  5 Discussion & Conclusion

To quantify the effect of possible missing faint galaxies on the total IGL, we
calculated and summed the IGL fluxes in each filter and for each replication.
Then, using the fact that galaxies with 24 $\sim$$<$ $m_{\rm AB}$ $\sim$$<$ 29
mag make up at most 10% of the IGL (Windhorst et al., 2022), we used the
equation below to calculate the percent change in the IGL resulting from the
$2\times,$ $3\times$, and $4\times$ replications:

$\Delta\text{IGL}_{TOT}(n)=\frac{F_{n}-F_{1}}{10F_{1}}$ (3)

where $F_{n}$ is the integral of the flux for the $n$th replication, and
$F_{1}$ is the integral of the flux from the original XDF data in that given
filter. $F_{1}$ is multiplied by 10 simply because it represents 10% of the
IGL on its own. The denominator therefore represents 100% of the IGL. The
results of this analysis show that, in all filters, doubling the number of
detectable objects with 24 $\sim$$<$ $m_{\rm AB}$ $\sim$$<$ 29 mag increases
the total IGL by 8-10%, tripling the number increases the IGL by 16-20%, and
quadrupling the number increases it by 24-30%. For the F606W filter, this
means that the missing IGL due to natural confusion increases to $30\%$ as the
galaxy number counts increase from $\sim 1.3\times 10^{6}$ to $\sim 5.3\times
10^{6}$ objects per square degree.

To address how many galaxies could be missing from the current XDF data, we
extrapolated the linear trend in relative source counts to a $1\times$ value
for the three ACS/WFC filters. This calculation yields estimated total
completeness percentages at faint magnitudes of 93% for F435W filter, 99% for
the F606W filter, and 94% for the F775W filter.

Based on this analysis, it does not appear likely that missing faint galaxies
explain the optical EBL discrepancy. The completeness limits of the WFC3/IR
images analyzed here (see the Appendix) were slightly lower than that of the
ACS/WFC images. This is likely due to the larger FWHM at the wavelengths of
the WFC3/IR filters. See SVO Filter Profile Service and Windhorst et al.
(2011, their Table 2) for the FWHMs and more info on these _HST_ filters. The
WFC3/IR filters were analyzed in this work mostly to ensure that the analysis
was consistent, and to determine there was no major wavelength or FWHM
dependence in our results. The quoted fractions in the legends of Figures 3
and 5 are the result of the total number of objects with 24 $\sim$$<$ $m_{\rm
AB}$ $\sim$$<$ 29 recovered in the denser realizations compared to the
original after appropriate normalization. One can still analyze completeness
at specific magnitudes based on these results. For example, in most filters
the downward trend in the ratio points does not begin until $m_{\rm AB}$
$\sim$ 27 mag, and most filters were at least 50% complete up to $m_{\rm
AB}$$=29$ mag, even in the $4\times$ realizations. In other HUDF works, object
counts are also done in other filters, but we confined ourselves to the six
deepest HUDF filters, which exclude F814W, F850LP, F140W, and the 3 UV filters
of Teplitz et al. (2013).

It is clear that significant numbers of faint galaxies are still missed or
blended with other objects in this experiment. This is expected because of the
large number of galaxies that are being added to the original, especially in
the 3$\times$ and 4$\times$ replications. Natural confusion and surface
brightness limitations cause the incompleteness of the sample to be noticeable
for sizes bigger than 0$\farcs$15–0$\farcs$2, and magnitudes $m_{\rm AB}$
$\sim$$>$ 29 mag. Even so, our conclusions are two-fold. If there existed 4
times more faint galaxies in our Universe than what we observe today, then:
(1) we would only miss at most 30% of them (i.e., the highest,
$4\times$-replicated relative number count deficit in Figure 5), and so galaxy
catalogs would contain many more galaxies than are present in the actual XDF
data; and (2) not nearly enough IGL would be missed due to object overlap in
this sample (at most 30% of the total) to account for the 100% IGL-EBL
discrepancy.

Further studies of truly hierarchical simulations could be done to find out
how many, and what kind, of galaxies are missed on a portion of the sky when
their total surface density changes. The SKYSURF project (Windhorst et al.,
2022; Carleton et al., 2022) is working to resolve the EBL discrepancy by
reprocessing _HST_ data to include and measure the sky surface brightness for
both EBL and diffuse Zodiacal light calculations.

## Acknowledgments

We thank the anonymous referee for insightful comments that helped improve
this paper. We also thank the entire SKYSURF team for their interest,
comments, and discussion about this work; especially, we thank Rosalia O’Brien
for her support and helpful feedback. Finally, we thank Alex Van Engelen and
his research group at ASU for their insight and perspective. This project is
based on observations made with the NASA/ESA _Hubble Space Telescope_ and
obtained from the Hubble Legacy Archive, which is a collaboration between the
Space Telescope Science Institute (STScI/NASA), the Space Telescope European
Coordinating Facility (ST-ECF/ESA), and the Canadian Astronomy Data Centre
(CADC/NRC/CSA). We acknowledge support for _HST_ program AR-15810 provided by
NASA through a grant from STScI, which is operated by the Association of
Universities for Research in Astronomy, Incorporated, under NASA contract
NAS5-26555.

## References

  * Abeysekara et al. (2019) Abeysekara, A. U., Archer, A., Benbow, W., et al. 2019, The Astrophysical Journal, 885, 150.  http://dx.doi.org/10.3847/1538-4357/ab4817
  * Acciari et al. (2019) Acciari, V. A., Ansoldi, S., Antonelli, L. A., et al. 2019, Monthly Notices of the Royal Astronomical Society, 486, 4233–4251.  http://dx.doi.org/10.1093/mnras/stz943
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, Nature, 440, 1018
  * Beckwith et al. (2006) Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729
  * Bernstein (2007) Bernstein, R. A. 2007, ApJ, 666, 663
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
  * Bunker et al. (2010) Bunker, A. J., Wilkins, S., Ellis, R. S., et al. 2010, Monthly Notices of the Royal Astronomical Society, 409, 855.  https://doi.org/10.1111/j.1365-2966.2010.17350.x
  * Carleton et al. (2022) Carleton, T., Windhorst, R. A., O’Brien, R., et al. 2022, SKYSURF: Constraints on Zodiacal Light and Extragalactic Background Light through Panchromatic HST All-Sky Surface-Brightness Measurements: II. First Limits on Diffuse Light at 1.25, 1.4, and 1.6 microns, arXiv, doi:10.48550/ARXIV.2205.06347.  https://arxiv.org/abs/2205.06347
  * Casertano et al. (2000) Casertano, S., de Mello, D., Dickinson, M., et al. 2000, AJ, 120, 2747
  * Dressel (2022) Dressel, L. 2022, Wide Field Camera 3 Instrument Handbook, Version 14.0, (Baltimore: STScI)
  * Driver (2021) Driver, S. P. 2021, arXiv e-prints, arXiv:2102.12089
  * Driver et al. (2016) Driver, S. P., Andrews, S. K., Davies, L. J., et al. 2016, The Astrophysical Journal, 827, 108.  https://doi.org/10.3847/0004-637x/827/2/108
  * Ellis et al. (2013) Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7
  * Finkelstein et al. (2012) Finkelstein, S. L., Papovich, C., Ryan, R. E., et al. 2012, The Astrophysical Journal, 758, 93.  https://doi.org/10.1088/0004-637x/758/2/93
  * Hill et al. (2018) Hill, R., Masui, K. W., & Scott, D. 2018, Applied Spectroscopy, 72, 663–688.  http://dx.doi.org/10.1177/0003702818767133
  * Illingworth et al. (2013) Illingworth, G. D., Magee, D., Oesch, P. A., et al. 2013, ApJS, 209, 6
  * Joye & Mandel (2003) Joye, W. A., & Mandel, E. 2003, in Astronomical Society of the Pacific Conference Series, Vol. 295, Astronomical Data Analysis Software and Systems XII, ed. H. E. Payne, R. I. Jedrzejewski, & R. N. Hook, 489
  * Keenan et al. (2010) Keenan, R. C., Barger, A. J., Cowie, L. L., & Wang, W.-H. 2010, The Astrophysical Journal, 723, 40–46.  http://dx.doi.org/10.1088/0004-637X/723/1/40
  * Koekemoer et al. (2013) Koekemoer, A. M., Ellis, R. S., McLure, R. J., et al. 2013, ApJS, 209, 3
  * Kron (1980) Kron, R. G. 1980, ApJS, 43, 305
  * Lauer et al. (2021) Lauer, T. R., Postman, M., Weaver, H. A., et al. 2021, The Astrophysical Journal, 906, 77.  http://dx.doi.org/10.3847/1538-4357/abc881
  * Lauer et al. (2022) Lauer, T. R., Postman, M., Spencer, J. R., et al. 2022, The Astrophysical Journal Letters, 927, L8.  https://doi.org/10.3847/2041-8213/ac573d
  * Madau & Pozzetti (2000) Madau, P., & Pozzetti, L. 2000, Monthly Notices of the Royal Astronomical Society, 312, L9–L15.  http://dx.doi.org/10.1046/j.1365-8711.2000.03268.x
  * Matsumoto et al. (2015) Matsumoto, T., Kim, M. G., Pyo, J., & Tsumura, K. 2015, The Astrophysical Journal, 807, 57.  https://doi.org/10.1088/0004-637x/807/1/57
  * Matsuoka et al. (2011) Matsuoka, Y., Ienaka, N., Kawara, K., & Oyabu, S. 2011, ApJ, 736, 119
  * Mattila & Väisänen (2019) Mattila, K., & Väisänen, P. 2019, Contemporary Physics, 60, 23
  * Oesch et al. (2009) Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2009, The Astrophysical Journal Letters, 0909.5183
  * Oesch et al. (2010) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2010, ApJ, 709, L16
  * Rafelski et al. (2015) Rafelski, M., Teplitz, H. I., Gardner, J. P., et al. 2015, arXiv, 1505.01160
  * Roseboom et al. (2010) Roseboom, I. G., Oliver, S. J., Kunz, M., et al. 2010, MNRAS, 409, 48
  * Ryon (2022) Ryon, J. E. 2022, ACS Data Handbook for Cycle 30, Version 21.0, (Baltimore: STScI)
  * Serjeant et al. (1997) Serjeant, S. B. G., Eaton, N., Oliver, S. J., et al. 1997, Monthly Notices of the Royal Astronomical Society, 289, 457.  https://doi.org/10.1093%2Fmnras%2F289.2.457
  * Silva et al. (2005) Silva, L., De Zotti, G., Granato, G. L., Maiolino, R., & Danese, L. 2005, Mon. Not. Roy. Astron. Soc., 357, 1295
  * Teplitz et al. (2013) Teplitz, H. I., Rafelski, M., Kurczynski, P., et al. 2013, The Astronomical Journal, doi:10.1088/0004-6256/146/6/159
  * Williams et al. (1996) Williams, R. E., Blacker, B., Dickinson, M., et al. 1996, AJ, 112, 1335
  * Windhorst et al. (2008) Windhorst, R., Hathi, N., Cohen, S., et al. 2008, Advances in Space Research, 41, 1965
  * Windhorst et al. (2011) Windhorst, R. A., Cohen, S. H., Hathi, N. P., et al. 2011, ApJS, 193, 27
  * Windhorst et al. (2021) Windhorst, R. A., Carleton, T., Cohen, S. H., et al. 2021, Galaxy Science with ORCAS: Faint Star-Forming Clumps to AB$\leq$31 mag and re$\geq$ 0.01”, arXiv, doi:10.48550/ARXIV.2106.02664.  https://arxiv.org/abs/2106.02664
  * Windhorst et al. (2022) Windhorst, R. A., Carleton, T., O’Brien, R., et al. 2022, SKYSURF: Constraints on Zodiacal Light and Extragalactic Background Light through Panchromatic HST All-Sky Surface-Brightness Measurements: I. Survey Overview and Methods, arXiv, doi:10.48550/ARXIV.2205.06214.  https://arxiv.org/abs/2205.06214

## Appendix

Below are the equivalent plots to those in Figure 3 for the remaining XDF
filters mentioned above. Trends are similar, though for the WFC3/IR filters
the total percent recovery of object counts are $\sim$10-20% lower than those
in the ACS/WFC filters due to the filters’ higher FWHMs, as discussed in § 5.

Figure 5: _Left Column:_ Histograms of object counts per $\frac{1}{2}$ mag bin
per deg2 for the F105W, F125W, and F160W original and replicated images up to
their completeness limits. AX Rep refers to the integer A replications of the
XDF data, and AxXDF refers to that integer times the original XDF bin counts.
See Figure 3 for the ACS/WFC filters. The values in the legend are up to the
completeness limit in each filter and replication. _Right Column:_ The ratios
per $\frac{1}{2}$ mag bin of the replicated image bin counts over two, three,
and four times the original image bin counts. These curves demonstrate the
completeness of each replication and how it varies with magnitude.