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

Extragalactic background light (EBL) is the dominant background in the universe after the cosmic microwave background, and is generally considered to consist of the short wavelength (0.1 μm) ultraviolet to long wavelength (1000 μm) far-infrared emission (Driver et al. 2016; Mattila & Vaisanen 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, active galactic nuclei, 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 galaxy 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 the 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-rays and the EBL photons interact, constraining the permitted levels of the EBL intensity with a ∼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 that 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_AB > 24 mag, i.e., that surveys miss significant numbers of faint galaxies with m_AB ≲ 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 ≲ m_AB ≲ 30 mag to make up for the missing diffuse EBL flux (Bernstein 2007; Hill et al. 2018; Lauer et al. 2021, 2022). Deep extragalactic surveys with the Hubble Space Telescope (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; Illingworth et al. 2013; Koekemoer et al. 2013; Rafelski et al. 2015) have detected and resolved distant galaxies, with some candidates above z ∼ 10 and to apparent magnitudes m_AB ∼ 30 mag, constraining the observed galaxy luminosity function (Oesch et al. 2009; Windhorst et al. 2011, 2021; Finkelstein et al. 2012). 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).

This is the third Letter associated with the SKYSURF project, an Archival HST program (AR-15810; PI: R. Windhorst) that aims to resolve the direct EBL versus IGL conflict. SKYSURF is addressing the conflict by reprocessing HST data in order to measure the sky values, instead of reducing them (Carleton et al. 2022; O'Brien et al. 2022; Windhorst et al. 2022)

Windhorst et al. (2022) showed that bright objects with m_AB < 17 mag contribute ∼25% of the IGL, while objects 17 ≲ m_AB ≲ 22 mag contribute about 50% and objects fainter than m_AB = 22 mag contribute 25% of the IGL. Moreover, objects fainter than m_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_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 Section 2 we discuss the HST data used in this work and how it was preprocessed. Section 3 discusses the use of SourceExtractor on these images. In Section 4 we show and analyze the SourceExtractor source photometry extracted from the deep field images, and in Section 5 the conclusions and implications of these results are laid out.

It was desired to perform this exercise with data rather than simulations (such as those in Windhorst et al. 2022 Section 4.2) so that galaxy sizes and morphologies could be as realistic as possible. We therefore use the XDF data set 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 ∼0.12–0.2 mag 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 broadband filters, but more filters were used in the analysis (F105W, F125W, and 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: 10.17909/T9RG6J. 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 × 1401 pixel region centered on pixel (2734,3072). Figure 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.

Zoom In Zoom Out Reset image size

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 (Figure 1(b)). Each excised ∼2870 × 2870 pixel region is centered on pixel (2625,2625) in the original XDF mosaics.

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_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_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 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× the sky rms value (which is ∼53% lower than the 1.5σ 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.

Zoom In Zoom Out Reset image size

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 coadded to the original mosaics. This minimizes the negative effects of coadding sky noise.

When the ACS/WFC and WFC3/IR images were drizzled from their native pixels to the 60 mas 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 × 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 × 90°. Each "rotated" image was coadded to the original. This multiplied the number of potentially detectable objects by 2×, 3×, or 4×, 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× replication, the weight image was constructed from the original XDF weight map and the rotated weight maps:

$$\begin{eqnarray}&&{\rm{W}}3{\rm{x}}=\displaystyle \frac{1}{1/{W}_{0}+1/{W}_{90}+1/{W}_{180}},\end{eqnarray} \tag{ 1 }$$

where W₉₀ 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σ,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 Figures 3, 4, and 5. 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 Figure 3), the error bars were calculated simply by ${\sigma }_{\mathrm{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 Figure 3), the uncertainties are based on their corresponding histograms and derived for a Gaussian error distribution:

$$\begin{eqnarray}&&{\sigma }_{\mathrm{bin}}=\displaystyle \frac{m}{{pn}}\ast \sqrt{{\left(\displaystyle \frac{\sqrt{m}}{m}\right)}^{2}+{\left(\displaystyle \frac{\sqrt{{pn}}}{{pn}}\right)}^{2}},\end{eqnarray} \tag{ 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 coadded images. We note that the rms noise increased by 4%–7% from the initial XDF mosaics to their 4× replications. This can be explained by the wings of some objects going undetected by SourceExtractor, and therefore being included in sky statistics.

Zoom In Zoom Out Reset image size

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 the 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).

We now wish to examine if the transposed and coadded images do indeed contain 2×, 3×, or 4× as many objects as the original, respectively. 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_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×, 3×, or 4× 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_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_AB ≃ 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.

Zoom In Zoom Out Reset image size

Comparing slopes in the left column of Figures 3–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 Section 5.

It is worth noting that the exercises in this Letter were initially performed on unmasked 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 unmasked replications showed that the unmasked analysis had ∼5%–30% lower total relative number counts (right column of Figures 3 and 5), meaning fewer objects were recovered in the unmasked 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 ≲ m_AB ≲ 29 mag were on average 10% lower than those from the unmasked images. This is due to the generally larger number of objects in the unmasked image. Therefore, our completeness curves and summed values in each HST filter, in the right panels of Figures 3 and 5, are upper limits.

Zoom In Zoom Out Reset image size

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 ≲ m_AB ≲ 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×, 3×, and 4× replications:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{IGL}}_{\mathrm{TOT}}(n)=\displaystyle \frac{{F}_{n}-{F}_{1}}{10{F}_{1}},\end{eqnarray} \tag{ 3 }$$

where F_n is the integral of the flux for the nth replication, and F₁ is the integral of the flux from the original XDF data in that given filter. F₁ 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 ≲ m_AB ≲ 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 ∼1.3 × 10⁶ to ∼5.3 × 10⁶ 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× value for the three ACS/WFC filters. This calculation yields the estimated total completeness percentages at faint magnitudes of 93% for the 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 Spanish Virtual Observatory 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 that 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 ≲ m_AB ≲ 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_AB ∼27 mag, and most filters were at least 50% complete up to m_AB = 29 mag, even in the 4 × 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 three 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× and 4× replications. Natural confusion and surface brightness limitations cause the incompleteness of the sample to be noticeable for sizes bigger than 015–02, and magnitudes m_AB ≳ 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× -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 (Carleton et al. 2022; O'Brien et al. 2022; Windhorst et al. 2022) is working to resolve the EBL discrepancy by reprocessing the HST data to include and measure the sky surface brightness for both the EBL and diffuse Zodiacal light calculations.

We thank the anonymous referee for insightful comments that helped improve this Letter. 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. RAW, SHC, and RAJ acknowledge support from NASA JWST Interdisciplinary Scientist grants NAG5-12460, NNX14AN10G, and 80NSSC18K0200 from GSFC.

Software: SourceExtractor (Bertin & Arnouts 1996).

Figure 5 contains 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 ∼10%–20% lower than those in the ACS/WFC filters due to the filters' higher FWHMs, as discussed in Section 5.