Unresolved z~8 point sources and their impact on the bright end of the galaxy luminosity function

The distribution and properties of the first galaxies and quasars are critical pieces of the puzzle in understanding galaxy evolution and cosmic reionization. Previous studies have often excluded unresolved sources as potential low redshift interlopers. We combine broadband color and photometric redshift analysis with morphological selections to identify a robust sample of candidates consistent with unresolved point sources at redshift $z\sim8$ using deep Hubble Space Telescope images. We also examine G141 grism spectroscopic data to identify and eliminate dwarf star contaminants. From these analyses, we identify three, bright ($M_{UV}\lesssim-22$ ABmag) dropout point sources at $7.5<z<8.1$. Spectral energy distribution analyses suggest that these sources are either quasars or compact star-forming galaxies. The flux captured by the IRAC 4.5 $\mu$m channel suggests that they have moderate $H\beta$+$[OIII]$ equivalent widths. We calculate the number density of point sources at $z\sim7$-8, and find that a double powerlaw model well describes the point source distribution. We then extend our analysis to estimate the combined point source + galaxy luminosity function and find that the point sources have a non-negligible contribution to the bright-end excess. The fact that the point sources dominate only at $M_{UV}\lesssim-22$ suggests that their contribution to cosmic reionization is likely limited. While spectroscopic follow-up is needed to confirm the nature of these point sources, this work demonstrates that the inclusion of Lyman dropout point sources is necessary for a complete census of the early galaxies at the epoch of cosmic reionization.


INTRODUCTION
Statistical studies of the first galaxies and quasars are crucial to understanding their formation and evolution processes. To-date, a tremendous amount of effort has been made to probe the early universe with high redshift surveys like CANDELS (Koekemoer et al. 2011;Grogin et al. 2011;Bouwens et al. 2019), BoRG (Trenti et al. 2011;Bradley et al. 2012;Morishita et al. 2018; Morishita 2021), HUDF12 , XDF (Illingworth et al. 2013), CLASH (Postman et al. 2012), HFF (Lotz et al. 2017), RELICS (Coe et al. 2019), ULTRA-VISTA (McCracken et al. 2012;Stefanon et al. 2017aStefanon et al. , 2019Bowler et al. 2020), among others. These surveys combined with follow-up spectroscopy have successfully identified some of the earliest galaxies up to z ∼ 9 -10, yet characterizing the number densities and physical properties of these early sources remain incomplete. It is necessary to accurately quantify the early populations with observational constraints.
Characterizing the luminosity function is a fundamental step in estimating the contribution from various luminous sources; it describes the number density of sources as a function of luminosity, or absolute magnitude. Since ultraviolet (UV) emission is primarily dominated by ionizing sources, the rest-frame UV luminosity function is a useful tool in investigating the early galaxy populations. In particular, the shape of the luminosity function can provide insights into the different physical processes such as star-formation and quasar activity that drive galaxy formation.
The faint end of the luminosity function is believed to be the key driver for cosmic reionization (e.g., Ishigaki et al. 2018;Atek et al. 2018), in which the early universe transitioned from completely neutral to almost ionized (Ouchi et al. 2010;Konno et al. 2014;Pentericci et al. 2014;Robertson et al. 2015;Mason et al. 2018Mason et al. , 2019Hoag et al. 2019). It is believed that reionization paved the way for the formation of the first galaxies (Loeb & Barkana 2001), yet the question of what astrophysical objects are primarily responsible for reionization remain debated.
The bright end of the luminosity function is composed of the brightest sources that may be signposts of in-situ star formation or even quasar activity. The discovery of luminous quasars (Mortlock et al. 2011;Bañados et al. 2018a;Yang et al. 2020;Wang et al. 2021a) and luminous star forming galaxies at z ∼ > 7 (Zitrin et al. 2015;Oesch et al. 2016;Hashimoto et al. 2018;Jiang et al. 2021) may indicative of populations of luminous sources that are unaccounted for. There is no consensus on the shape of the bright end luminosity function; some even suggest that the early galaxy luminosity function may depart from the standard Schechter form (Harikane et al. 2022). This departure manifest as a bright end excess (e.g., Morishita et al. 2018), and its origins remain unclear. Theoretical studies suggest that this bright excess may be caused by intense and compact star-forming clumps (e.g., Ma et al. 2018) or even stochastic quasar activity (e.g., Ren et al. 2020). Luminous source are also believed to contribute to cosmic reionization to some degree, yet the consensus on their contribution remain controversial (e.g., Willott et al. 2010;Finkelstein et al. 2015;Jiang et al. 2016;Matsuoka et al. 2019;Naidu et al. 2020). The characterization of the brightest sources at high redshifts remain elusive.
High-redshift sources are typically identified with the Lyman dropout photometric selection (Steidel et al. 1996) combined with follow-up spectroscopic confirmation. However, the complication is that the redshifted spectral energy distribution (SED) of these sources at the end of cosmic reionization, z ∼ 7-8, overlaps with those of low-mass foreground stars. As a result, previous studies have often excluded compact, unresolved sources with star-like morphology in preference for more galaxy-like sources with extended morphology. Nevertheless, there is evidence from lensing surveys that early galaxies are very compact (e.g., Bouwens et al. 2017;Salmon et al. 2020). Some even predict compact starforming clumps (Ma et al. 2018). So there is a possibility that certain population of quasars and compact galaxies at at z ∼ 7-8 that are rejected with the standard selection.
There is a renewed interest in examining these overlooked point sources. A recent medium-depth, wide ∼ 0.4 deg 2 HST survey (SuperBoRG; Morishita et al. 2020;Morishita 2021) has identified several z 8 point sources as potential quasar candidates. Key features of these z ∼ 8 point sources are their blue rest-frame UV slopes and the Spitzer /IRAC flux excess in the SED, which may be indicative of significant Hβ and [O iii] emission often seen in quasars. Observations suggest that these point sources are unlikely to be foreground stars and may contribute to the bright end of the luminosity function.
In this study, as part of our HST archival program (AR 15804; PI. Morishita), we reexamine the selection of high-redshift compact unresolved (point) sources that have been overlooked in previous studies. We take advantage of the successful SuperBoRG study to revisit z ∼ 7-8 point sources in the CANDELS legacy HST fields. Since the selection criteria of our study and of SuperBoRG are complementary, we combine the results of both studies to characterize the z ∼ 7-8 dropout point sources and to quantify their contribution to the total galaxy luminosity function. For simplicity, we will refer to these sources as "point sources" throughout the paper.
This paper is organized as follows. In Section 2 we describe the data reduction and target selection from 3D-HST. In Section 3 we explore the properties of the targets selected. And in Section 4 we discuss the physical implications of these objects. We use the AB-magnitude system (Oke & Gunn 1983;Fukugita et al. 1996) and adopt the h = 0.7, Ω M = 0.3, and Ω Λ = 0.7 cosmology.

TARGET SELECTION 2.1. Source Catalog
Our primary focus is to identify sources that satisfy the dropout color selection and have point source morphology in the CANDELS fields. We begin our analyses with the publicly available photometric catalogs provided by the 3D-HST team (Brammer et al. 2012;van Dokkum et al. 2013). The 3D-HST is a HST nearinfrared spectroscopic survey designed to study galaxies across the universe. It surveyed nearly 700 arcmin 2 of the well-studied HST /CANDELS Treasury fields to obtain direct images and spectroscopic data with the ACS/G800L and WFC3/G141 grisms. 3D-HST covers about 75% of the original CANDELS area. When all the fields are combined, their photometric observations of H 160 reach median 5 σ depths at 26 mags at 1 aper-ture. Further details of the survey and the published catalog can be found in Skelton et al. (2014); Momcheva et al. (2016).
Our choice of using the 3D-HST catalog over the catalogs published by the CANDELS team (Guo et al. 2013;Galametz et al. 2013;Stefanon et al. 2017b;Nayyeri et al. 2017;Barro et al. 2019) is that uniform analysis is performed on all 5 CANDELS fields by the 3D-HST team to create the source catalog. This vastly simplifies the source detection procedure (Sec. 2.2) and the completeness simulation analysis (Sec. 3.3), to calculate the number density of the target population. However, due to inconsistencies in the filter coverage, we only analyze 4 of the 5 CANDELS fields (AEGIS, COSMOS, GOODS South, and UKIDSS-UDS), where F814W, F125W, and F160W filters are available (Sec. 2.2). We also exploit the published G141 grism data to identify low redshift interlopers (Sec. 2.4).
We obtain deep HST data from the publicly available 3D-HST database. The HST image mosaics used have already been corrected for distortions and drizzled to the plate scale of 0. 06 pixel −1 . The photometric source catalogs were produced using pointspread function (PSF)-matched aperture photometry, reduced using SExtractor (Bertin & Arnouts 1996), and flux calibrated to an aperture radius of 0. 7. The HST ACS and WFC3 images were convolved to match the HST /F160W PSF (∼ 0. 14). Ground-based optical, NIR, and Spitzer /IRAC fluxes are similarly PSFmatched to a combination of F125W, F140W, and F160W priors and aperture corrected to F160W (or F140W, otherwise).

Color-dropout and shape selection
Our strategy in identifying z ∼ 7-8 point source candidates from the 3D-HST is twofold. First, we identify sources with the Lyman-break dropout technique (Steidel et al. 1996) from the photometric catalog. Then, we select point sources from the list of Lyman-dropout sources. The color selection is only based on deep HST photometry. A caveat is that unlike the Morishita et al. (2020); Morishita (2021) selection, which uses the F105W/F125W/F160W (Y 105 /J 125 /H 160 ) filters, the 3D-HST catalog does not include Y 105 fluxes. Instead, we use the  color-dropout criteria, which is based on the F814W/F125W/F160W (I 814 /J 125 /H 160 ) selection: Compared to the Morishita et al. (2020) selection, this color selection results in a broader z ∼ 7-8 selection.  Table 1 are shown as red stars. Point sources that do not meet the photometric redshift selection are shown as dark red stars. Spectroscopically confirmed dwarf stars are shown as gray diamonds. We also compare with dwarf star template colors (Burgasser 2014) in gray dots, predicted quasar colors from composite SDSS quasar spectrum (Pâris et al. 2018) in blue, and simple powerlaw spectra at β = −1, −2, −3 in green. Strong overlap of known dwarf star colors and our Lyman-dropout sources indicate the need for further follow-up to disentangle the degeneracy.
Also, the GOODS-North catalog does not include I 814 data, so we only examine 4 of the 5 CANDELS fields. Where available, we use the HST /ACS blue filters (I 814 and bluer) to determine the Lyman-dropout with strict blue signal-to-noise constraints, combined with the nondetection flag from the HST pipeline catalog. Since we require 2 σ non-detections for I 814 fluxes, we calculate the resulting I 814 − J 125 lower limit color as: (I 814 − J 125 ) lim = −2.5 log 10 (2σ 814 /J 125 ). (2) We do not include sources with I 814 −J 125 that fall below the non-detection limit (i.e. falls outside the selection box in Figure 1). While additional ground-based fluxes are available, higher S/N HST fluxes are prioritized for color selection here (but see Sec. 2.4). From the color dropouts, we identify point sources based on two morphological parameters, elongation and flux concentration, measured in the H 160 filter. The selection criteria are defined by Morishita et al. (2020): Table 1. Target selected from the 3D-HST fields based on color and morphology criteria. Photometric redshifts estimated with EAZY. The targets listed here are predicted to have photometric redshift probability of p(z ph > 6) > 70%. We also indicate whether grism spectra is available in the 3D-HST database. Fluxes are listed in Table 2. We also show the EAZY χν goodness-of-fit for the z ph and dwarf star SED fits.  e < 1.2 f 5 /f 10 > 0.5. ( Elongation, e (ratio of semi-major/semi-minor axes), describes the circularity of the source, and the flux concentration (flux ratio between inner and outer radii) describes the compactness. Morishita (2021) find light concentration is an appropriate metric for point source selection. We obtain e from the 3D-HST catalog. To calculate flux concentration, we run SExtractor on the image mosaics, matching the 3D-HST detection parameters, and obtain detailed aperture photometry of the targets. We extract the aperture fluxes to calculate the flux ratios. After careful comparison of the different H 160 flux ratios at different radii, we determine the f 5 /f 10 flux concentration, the flux ratios taken within the 5 pixel (0. 3) and 10 pixel (0. 6) radii, to be the appropriate criteria. This decision is based on the ability to concurrently recover known dwarf star contaminants due to the point source selection (see Sec. 2.4). Although the 3D-HST source catalogs include the star class flag that classify whether a source is star-like, Finkelstein et al. (2015) and Morishita (2021) have demonstrated that these flags are not complete down to fainter magnitudes; it fails to distinguish between fuzzy circular objects and compact point sources in the faint magnitude ranges up to ∼ 24 mag. We note that other studies (e.g., Roberts-Borsani et al. 2016; have successfully identified sources by combining color dropout, stellarity parameters, and SED properties. Our aim is to identify additional sources that may be missed with the standard method.
Of the 169,614 objects listed in the 3D-HST catalog, we identify 22 I 814 /J 125 /H 160 dropout point sources. Of these point sources, 7 meet the photometric redshift selection of z ph > 7, discussed in Section 2.3. We show these J 125 − H 160 vs. I 814 − J 125 color-color diagram of all color and z ph > 7 selected point sources in Figure  1. Then, we check the grism spectra to further eliminate any low redshift interlopers, discussed in Section 2.4. Our final z ∼ 7-8 point source candidates list consists of 3 point sources, which are listed in Table 1. We also visually inspect the HST images of the z ph selected Note- † The Y105 flux for EGS 29337 is taken from the Y105 − J125 value in Roberts-Borsani et al. (2016). ‡ Spitzer /IRAC fluxes with visually confirmed contamination are excluded. Table 3. Bagpipes SED fit priors assuming a young stellar population model. The model redshifts are fixed to the EAZYderived z ph found in Table 1. τage is the range of universe ages calculated with the exponential star-formation history model; M is the final stellar mass formed, Z is the metallicity, AV measures the dust attenuation, and U is the nebular ionization parameter that captures the nebular emission and continuum components.

Parameter Units Priors
targets to eliminate image artifacts and/or other spurious detections; postage stamp images of the final sample are shown in Figure 2 and discussed in Section 4. The observed fluxes of the point sources are shown in Table 2. We also list low confidence, I 814 − J 125 limited, non-candidate point sources in Appendix A.

Photometric redshifts and SED fits
The Lyman-break color selection is comprehensive but also allows low-redshift sources with similar colors to migrate into the selection window. To filter out these contaminants, we apply further selection based on the photometric redshift measurement discussed here.
Following the similar approach by Roberts-Borsani et al. (2021), we estimate the photometric redshift, z ph , using the photometric redshift code, EAZY (Brammer et al. 2008). While photometric redshifts are also included in the public 3D-HST catalog, they were calculated with a maximum limiting redshift of z = 6. Hence, we re-reduce the redshift estimates for all of our color-selected samples. We use EAZY in the default setup (v1.3 templates) to derive the best-fit SED and redshift posterior probability, p(z ph ). To ensure a more accurate redshift derivation, we use all available photometric data points, including ground-based fluxes that were excluded in our initial color selection in Sec. 2.2. We turn off magnitude priors in the fit to avoid any biased redshift selections, and fit between 0 ≤ z ph ≤ 9. The z ph fit range is based on the redshift probability from the survey completeness (Sec. 3.3). From this analysis, we eliminate low redshift interlopers by making a cut in which the probability of z ph > 6 is greater than 70%. The z ph and corresponding probability of our targets are shown in Table 1.
Upon determining z ph , we refine the SED fits using Bagpipes (Carnall et al. 2018) to determine their physical properties. Morishita et al. (2020) notes that distinguishing between luminous galaxies and quasars at z ∼ 7-8 is ambiguous and challenging without spectroscopy. In Figure 3 we plot both the best-fit quasar (described in Appendix B) and star-forming galaxy SEDs (from Bagpipes fitting described next), which clearly show degenerate profiles, with the exception of EGS 29337 (to be discussed later). Since precise modeling is beyond the scope of this study, in this paper we instead assume that the sources are well represented with a young stellar spectrum with nebular emission with Bagpipes modeling and explore the inferred properties.
We describe our Bagpipes fit methodology here. The redshift is fixed to z ph from EAZY, and we freely fit for other properties. The model fit priors are listed in Table 3, following the treatment in Roberts-Borsani et al.
(2021) as a guide. The best-fit Bagpipes SEDs and EAZY z ph probability distributions are shown in Figure  3. For each source, we use the best fitting SED model to extract the source's rest-frame UV luminosity, M U V , and other stellar parameters of interest for subsequent analysis, which is discussed in Section 3. The best-fit Bagpipes model parameters are shown in Table 4.  Table 3. PSF-limited R eff are shown as upper limits.

Excluding low-redshift contaminants
Figure 3. Best-fit Bagpipes galaxy SEDs (black) and the 1σ distribution (light gray) based on the broadband photometric data (red diamonds) with 1σ uncertainties. The x-axis errorbars indicate the filter bandwidth. We note that the we enforce 5σ cuts on HST J125 and H160 fluxes. The large J-band fluxes, relative to J125 fluxes, in EGS 515 and EGS 29337 correspond to shallower ground-based CFHT data (Bielby et al. 2012). The z ph are fixed based on the best EAZY estimates. The non-detections are indicated with 2σ upper limits. We plot the best-fit quasar template SEDs (light blue). Quasars and starbursts are nearly degenerate at the predicted z ph . We can see that a galaxy model is most appropriate for EGS 29337, as expected from Roberts-Borsani et al. (2016); Stark et al. (2017). We also plot the best-fit dwarf star templates in dark gray, which also show nearly degenerate fits. The insert plot shows the EAZY z ph probability distribution in blue. The resulting p(z ph ) distribution across 6.5 z 8.5 is broad since Y105 data was not included.
To further exclude low redshift interlopers among the selected point sources, we utilize G141 grism spectra made available by the 3D-HST team. As alluded to earlier, dwarf star SEDs have a sharp 1 µm drop-off that mocks the Lyman break of z ∼ 7-8 objects, making them likely interloper contaminants. There are notable spectral features at 1 µm, 1.25 µm, and 1.6 µm that are captured by the G141 grism. We find that some of our point sources are not listed in the grism catalog, due either to the extraction limit (JH 140 ≈ 26 mag) or incomplete spectral coverage (landing on/outside of detector edge). It is noted that the G141 grism does not cover the redshifted Lyα break at z ∼ 7-8, making the confirmation as high redshift sources difficult. Therefore, the objective of our inspection here is to exclude interlopers through the detection of continuum spectral features instead of characterizing their spectra.
When extracted grism spectra are available for the sources selected in Sec. 2.2, we perform spectral fits to low-mass L and T dwarf template spectra, which were observed with the SpeX spectrograph on NASA InfraRed Telescope Facility (Rayner et al. 2003). Template spectra are obtained from the SpeX Prism Library (Burgasser 2014). Based the spectral fits, we identified 4 T dwarfs with clear spectral features. Their 1D spectra are shown in Appendix C. In fact, 3 of these dwarf stars were also identified in a recent 3D-HST dwarf star study (Aganze et al. 2022). We exclude these targets from the final point source list. We note that none of our final redshift selected, point source candidates were identified by Aganze et al. (2022). However, the Aganze et al. (2022) selection was limited to spectra with S/N > 10, and thus their selected sources are all brighter than our final targets (H 160 ∼ < 24 mag). This may simply reflect the limitations of the 3D-HST grism data instead of differences in selection. We also fit the photometric SEDs with SpeX dwarf star templates using EAZY. In Table 1 we list the χ ν,BD , and in Figure 3 we show the best-fit models. When compared to the χ ν of z ph fits, we find that our final point sources are better constrained as z ∼ 7-8 sources.

Visual inspection of point sources
As the final step of our sample selection, here we examine the images to identify and to eliminate any spurious fluxes in HST and Spitzer. Once we eliminate false detections, we repeat and refine the EAZY z ph estimates and Bagpipes SED fits.
The postage stamp images of our final point sources in Table 2. Images are extracted for available deep HST (Grogin et al. 2011;Koekemoer et al. 2011;Skelton et al. 2014) and Spitzer (Dickinson et al. 2003;Ashby et al. 2013) that make up the 3D-HST catalogs. From the images, we clearly see the blue color dropouts, which are also reflected in the SED fits. Some sources show suspicious blue detections. For example, the GDS 29369 images show suspicious I 814 fluxes, despite meeting the non-detection criteria. After comparing the different filter images, we concluded that this is likely noise artifacts because their flux centroids do not match and the size is on the same order as the surrounding noise structure. The catalog also suggested spurious ground-based blue fluxes observed with Subaru (Taniguchi et al. 2007) for GDS 45797. However, careful examination of the Subaru images suggested that they are artifacts due to diffraction spikes from a nearby star. So, we treat these blue fluxes as non-detections in our analysis. Fortunately, the inclusion of these blue fluxes did not have a major affect on the z ph or the SED fits results. Another uncertainty comes from Spitzer /IRAC fluxes, which suffer from lower spatial resolution. The 3D-HST catalog includes IRAC contamination flags for each channel; however, the flags are based on contamination within 3 apertures, whereas our fluxes are taken at 0. 7 apertures. As a result, we individually inspected the IRAC images to decide whether or not to include them in our analysis. If we visually confirmed obvious contamination within a channel, we omitted its flux.
2.6. Cross-matching with Chandra X-ray catalog Finally, we also cross-matched our HST selected source with the deep X-ray Chandra catalogs for GOODS-South (Luo et al. 2017) and AEGIS (Nandra et al. 2015). Significant X-ray emission would strengthen the case for quasar candidacy. Previous tar- Table 5. Upper limits on the Chandra X-ray rest-frame fluxes and luminosities. We assume a simple powerlaw with Γ = 2 at z ph without any obscuration.
3. RESULTS: NATURE OF POINT SOURCES 3.1. On the point source selection We cross-match our point source selection with findings from , which examined the CANDELS, HUDF09, HUDF12, ERS, and BoRG/HIPPIES fields to estimate the galaxy UV luminosity function. Both studies appply similar Lymandropout color selections. The main difference is that we explicitly search for point-sources with the morphology selection defined in Sec. 2.2 using the f 5 /f 10 flux ratios, while  implements the stellarity parameter (e.g., star class flag), combined with SED fit photometry, to eliminate point-sources.
We find that only EGS 29337 is detected in both catalogs (EGSY-0120800269) with the separation of δr = 0. 11, which is well within the PSF uncertainty. In fact, this object has also been spectroscopically confirmed as as a z ph = 7.4 galaxy (Roberts-Borsani et al. 2016;Stark et al. 2017). This means that some near pointlike sources are in fact bona-fide galaxies. We discuss the implications in Section 4. Our SED modeling in Figure 3 also supports these observations. We note that the other sources were not identified by other z ∼7-8 studies Roberts-Borsani et al. 2016;Stark et al. 2017, e.g.). In fact, these sources have larger f 5 /f 10 values compared EGS 29337, which suggest that they appear more point-like, and thus more likely to have been rejected in the galaxy selection. Thus, a unique aspect of this study is that we explore the z ∼ 7-8 dropout point-sources that are often excluded in earlier studies of high-redshift galaxies.  Table 4.

Stellar SED fit properties of point sources
We list the best-fit stellar population properties of the point sources, which were estimated with Bagpipes SED fits, in Table 4. Due to the limitation of data, we assume that our detected point sources are well represented by a young stellar spectrum. The fit results predict sub-solar metallicities for nearly all of our point sources. Considering the young age of the universe , the metallicity estimates are not surprising. However, since uncertainties in broadband SED fits are strongly influenced by assumptions in the star formation history and age-metallicitydust degeneracies (e.g., Fig. 12 in Morishita et al. 2019), it is difficult to place any confident constraints on the metallicity evolution in the early universe. Instead we examine the size and star formation rate (SFR) density properties of the point sources.
We derive the projected physical size, R eff , from the half-light radius, which is calculated with SExtractor. When compared to the 3D-HST PSF limit, 0. 14, we find that 2 of the calculated R eff are upper limits. This is similar to the Morishita et al. (2020) results, as shown in Figure 4. This means that these sources are likely unresolved by HST, despite meeting our point source selection criteria.
Using the SFR estimated from the Bagpipes modeling, we also calculate the SFR density, Σ SFR , which is defined as the average SFR within a circle with radius R eff : The calculated Σ SFR serves as lower limits since its uncertainty is dependent on the upper limit uncertainty in R eff . We compare the inferred Σ SFR , which is calculated from the M U V -SFR relation (Kennicutt 1998;Ono et al. 2013) defined as follows: M U V = −2.5 log 10 Σ SFR · πR 2 eff 2.8 × 10 −28 (M yr −1 ) + 51.59.
(5) It appears that our point source candidates are highly compact star-forming objects. In Figure 4, we plot R eff and Σ SFR against their corresponding M U V . We also compare the Σ SFR redshift evolution. Our results appear to be consistent with trends of high redshift galaxies discussed in Ono et al. (2013); Holwerda et al. (2018), which predicts a greater number of compact galaxies with high star-formation at earlier epochs.
Using the best-fit Bagpipes SEDs, we calculate the UV continuum slope β U V . We adopt the formula defined (e.g., Dunlop et al. 2013): where J 125 and H 160 here are the best-fit magnitudes from the best-fit Bagpipes SEDs. The calculated β U V are listed in Table 4 with the mean slope ofβ U V = −1.90 ± 0.35. The resultingβ U V is consistent with β U V of known bright galaxies at z ∼7-8 (e.g., Dunlop et al. 2013;Bouwens et al. 2014). Lastly, we estimate the rest-frame equivalent width due to the Hβ + [O iii] emission lines from the bestfit SED for objects with sufficient IRAC fluxes. This is possible because Hβ + [O iii] emission from z ∼ 7-9 sources is well sampled by IRAC CH1 and CH2, 3.6 µm and 4.5 µm, (Roberts-Borsani et al. 2016). We calculate the equivalent widths as, where f cont is the underlying continuum flux obtained from the best-fit Bagpipes spectrum, f ch2 is the observed Spitzer /IRAC CH2 flux, ∆λ ch2 ∼ 1 µm is the full-width half maximum of the CH2 filter, and z ph is from EAZY. Comparing the values in Table 4 and the SED plots in Figure 3,

Number density of point sources
From the survey data, we constrain the point source luminosity function at z ∼ 7-8. We produce our own completeness simulation to calculate the effective volume, V eff , probed by the 3D-HST. We follow the completeness simulation treatment in Leethochawalit et al. (2021) to calculate V eff Carrasco et al. 2018;Calvi et al. 2016;Morishita et al. 2018).
We inject 500 sources into each 3D-HST image at each (M U V , z) bin: 100 ∆M U V bins across −26 ≤ M U V ≤ −16 and 13 ∆z bins across 7 ≤ z ≤ 9.4. All simulated sources in a given (M U V , z) grid have the same UV slope, which are randomly drawn from a Gaussian distribution with a mean slope ofβ U V = −2.2 ± 0.4. Source fluxes are calculated in the same way as in Skelton et al. (2014). We extract the simulated point sources according to our selection criteria described in Section 2.1. We repeat this process for every field. We show the redshift and magnitude probability distribution function of the extraction completeness in Figure 5.
After we calculate V eff , we estimate the number density of the point sources shown in Table 6 within each ∆M U V = 0.5 mag bin. We quote Poisson uncertainties for the number density (Gehrels 1986). In Figure 6, we plot the estimated number density of our point sources and compare with results from previous surveys of point sources and galaxies at z ∼ 7-8.
We fit the point source number density with both the Schechter function (Eq.8; Schechter 1976) and with the Double powerlaw (Eq.9; Hopkins et al. 2007), where φ * is the characteristic normalization, M * U V is the characteristic UV luminosity defined at M U V (1450Å), α defines the faint end slope, and β defines the bright end slope.
To properly include the bins of non-detection of the point source population at the faintest and brightest ends in fitting evaluation, we incorporate the upper limits from non-detections following the derivations from Sawicki (2012) for χ 2 minimization. The derivations for M U V and the χ 2 minimization are shown in Appendix Figure 5. The MUV and redshift dependent selection probability distribution as determined from our completeness simulation. The plot shows the ratio of color-selected point sources recovered to all input sources as a function of redshift and MUV . The colorbar to the right indicates this recovery fraction. At brighter MUV −21 magnitudes, at least 50% of simulated color-selected point sources are recovered (blue contour line). At fainter magnitudes, the recovery fraction decreases due to the combined effect of color-selection and source detection (gray contour lines at 40 % and 75 % detection). We also indicate the observed point source candidates at their respective z ph and MUV (red stars). D. We perform Markov Chain Monte Carlo (MCMC) sampling, using emcee package (Foreman-Mackey et al. 2013), to constrain the luminosity function. First we fit the luminosity function for point sources from this study (3D-HST selections). Then we apply the fits on all point sources selected from both this study and SuperBoRG.
If we freely fit for all of the parameters, the fit parameters do not converge to physically meaningful values. This is likely due to the lack of data points at both the faint and bright magnitude ranges. Instead, we opt for a more conservative approach and follow the known luminosity function shapes of z ∼ 7-8 galaxies and z ∼ 6 quasars. We fix both the faint end and bright end slopes and freely fit for φ * and M * U V . For the Schechter model fits, we fix the faint end slope to α = −2.2, which is the observed galaxy luminosity function at z ∼ 7-8 (Bouwens et al. 2021). For the double powerlaw mode fits, we fix both the faint end slope to α = −1.2 and the bright end slope to β = −2.7 based on the extrapolations of the quasar luminosity at z ∼ 6 (Matsuoka et al. 2018;Harikane et al. 2022).
With the deeper exposures of the 3D-HST survey, we improve the point source luminosity function fits estimated by Morishita et al. (2020). The best-fit luminos- Table 6. Number density, Φ, of 3D-HST point-sources at z ∼ 8. MUV is the rest-frame UV luminosity at 1450Å obtained from the best-fit Bagpipes SED. The V eff is calculated by the completeness simulation ). The 1σ uncertainties in Φ is based on Gehrels (1986). We also show the 1σ upper limits for Φ, where appropriate. We take ∆MUV ≈ 0.5 mag bins to match Bouwens et al. (2021).  Table 7. If we freely fit for M U V , we find that only the double powerlaw fit produces more reasonable parameters (M * U V ≈ −24) than the best-fit Schechter function, which instead suggests an unrealistically bright UV cutoff (M * U V ≈ −38, not shown). This may suggest that the point source luminosity function is more consistent with the high redshift quasar luminosity function (at z ∼ 6; Matsuoka et al. 2018). On the other hand, if we force a lower-bound on the M U V fit, then it is difficult to confidently favor either functions over the other. For both cases, the best-fit normalization φ * deviates from the Matsuoka et al. (2018) extrapolation by nearly ×100. This may be because 3D-HST is volume-limited, similar to the results in Morishita et al. (2020).

Number density of point source+galaxy populations
Finally, we fit the combined point source and galaxy luminosity functions. Here, we remove EGS 29337 from Table 7. Best fit luminosity function parameters. For the point sources, we fix both the α and β slopes depending on the model used and freely fit for φ * and M * U V . Parameters that are fixed are shown in a square bracket. The Schechter model slope is fixed to the galaxy faint-end slope from Bouwens et al. (2021), and the double power-law slopes are fixed to the AGN function slopes based on z ∼ 6 (Matsuoka et al. 2018;Harikane et al. 2022). We also fit for a combined galaxy and point source model with the slopes fixed and with freely fit parameters. The fit contours corresponding to the errors are shown in Figure 8.  Schechter or DPL). The main difference between the two studies is that Harikane et al. (2022) extrapolates the z ∼ 6 Matsuoka et al. (2018) quasar luminosity function to estimate the z ∼ 7 relation, while we use observed sources to anchor the point source luminosity function. Despite the differences in modeling, both results demonstrate that the inclusion of bright M U V −24 Super-BoRG sources at z ∼ 7-8 results in a departure from the previously measured galaxy luminosity function, i.e. a bright excess. We plot the z ∼ 7-8 point source number density and the best-fit luminosity functions in Figure  7. The best-fit parameters of all fits are shown in Table  7. The contours of the combined luminosity function fits are shown in Figure 8.
We compare the point source number densities and luminosity function against galaxy luminosity function measured at z ∼ 7-8 (Bouwens et al. 2021). We quantify the fraction of extended galaxies to all sources (galaxies and point sources) as a function of M U V as the following: where φ galaxy is the galaxy number density from the luminosity function (Bouwens et al. 2021) and φ points is the observed point source number density listed in Table  6. The uncertainty in f galaxy simply reflects the uncertainty in the number count of point sources (i.e. poisson). We plot f galaxy alongside the luminosity function fits in Figure 7 and compare with results from Harikane et al. (2022). We find that these point sources dominate at the bright M U V magnitudes. This suggests that the bright end excess implicated by the new point sources is not likely dominated by a typical population that has been identified in previous studies of high-redshift galaxies. We discuss the physical interpretations of this measured excess in the following section.
4. DISCUSSION Although our sources are selected with slightly different colors and from different surveys, the inferred number density of our point sources at M U V < −21.5 mags agree with the SuperBoRG point source study. This indicates that these z ∼ 7-8 dropout point sources are abundant enough to be detected in both surveys and are representative of similar populations. The inferred M U V suggests that these objects are driven by intense phenomena that occurs in a small physical scale, such as central starburst or quasar activity, which also shape their observed morphology point-like. In this section we explore the physical properties and implications of these sources.

Point sources as compact starburst galaxies
With the exception of EGS 29337, it is currently difficult to distinguish our final candidates between nonactive galaxies or quasar-hosting galaxies, the inferred sizes of the point sources is consistent with the observed trends of smaller galaxy sizes in the early universe (Oesch et al. 2010;Ono et al. 2013;Holwerda et al. 2015). If these point sources are compact non-active galaxies, we predict high Σ SF R based on our SED fitting analysis. Previously, Oesch et al. (2010); Ono et al. (2013) observed constant Σ SF R from z ∼ 4 to z ∼ 7 with a weak increase towards higher redshifts. Our SED analysis of the point sources appears to support this increasing Σ SF R trend, albeit we predict even larger values, as shown in Figure 4. Finally, given their predicted SFR and stellar masses, these sources may be progenitors of massive quiescent galaxies that are already present in the early universe (e.g., van Dokkum et al. 2008;Damjanov et al. 2009). This may suggest that our sources are UV enhanced starbursts and/or that additional physics may be at play. In fact, EGS 29337 has been shown to be one of the brightest z > 7 known with large SFR (Stark et al. 2017;Roberts-Borsani et al. 2016).

Point sources as low-luminosity quasars
While much of our SED analyses assumed the starforming SED properties, Figure 3 clearly shows that the SEDs of quasars and star-forming galaxies are degenerate. Theoretical predictions of early quasar properties suggest that variations in the quasar duty cycle may lead to an enhancement of UV bright quasars (Ren & Trenti 2021). With the detection of potentially UV bright quasars, there are also implications on the obscured quasar fraction, which is unknown at these redshifts (Vito et al. 2018(Vito et al. , 2019Inayoshi et al. 2020). Either there is an enhancement in the population of unobscured quasars or some physical mechanisms, such as  Table 7. The 1D histogram of the fit parameters are shown with the 16%, 50%, and 84% percentile values indicated.
powerful outflows, may drive obscured quasars to appear more luminous like unobscured quasars. Also, while non-detections from deep Chandra images suggest that no luminous quasar is present, the possibility of heavy Compton-thick obscuration may complicate this result (Ni et al. 2020). Another possibility is that these sources are quasars embedded in star-forming galaxies, similar to z ∼ 7 sources identified by Laporte et al. (2017). In fact, the possibility of either a compact starburst or quasar is supported by the recent discovery of an UV compact, red bright, X-ray faint object at z ∼ 7, which is hypothesized to be either a compact dusty, star-forming region or a Compton-thick super-Eddington quasar (Fujimoto et al. 2022).
Although we cannot distinguish between these possibilities with the current HST data, the James Webb Space Telescope's spatial resolution and sensitivity is expected reach below the predicted sizes and magnitudes of our sources at the range of −22 Marshall et al. (2021) predicts that deep NIRCam observations may allow us to study the quasar and its host galaxy at z ∼ 7. Its imaging and spectroscopic capability may enable us to confidently distinguish them as quasars or as compact star-forming galaxies. If they are revealed as quasars, they will become one of the most distant quasars ever discovered.

Investigating the impact on the bright-end excess
In this section, we focus on understanding the point sources' contribution to the bright end of the galaxy luminosity function. In Figure 7, we show the combined point source and galaxy luminosity function. Once the point sources from 3D-HST and SuperBoRG are incorporated, we find that the best-fit luminosity functions suggest the existence of a bright point source population, that may be missed by the galaxy surveys. This may support the existence of a bright end excess in the early universe. While ground-based observations (Bowler et al. 2020, e.g.) may detect unresolved sources, these studies are limited to select fields. Thus, we stress the importance of large-volume studies to accurately quantify the luminosity function.
Indeed, Harikane et al. 2022 present a comprehensive analysis of the luminosity function by combining quasar and galaxy populations identified in the HSC program. They proposed several explanations for the apparent bright end excess seen in galaxy luminosity functions. Physical mechanisms such as inefficient mass quenching due to high star-formation and/or poor quasar feedback, low dust obscuration in the host galaxy (Marques-Chaves et al. 2020, 2021, or even additional hidden quasar activity may increase the observed rest-frame UV luminosity (Mirocha 2020). Variations in the quasar duty cycle can also enhance the UV luminosity and contribute to the bright end (Ren & Trenti 2021). Our predictions as compact star-forming galaxies or quasars are consistent with these possibilities.
Other possibilities include superposition of lensed galaxies or even merging galaxies. However, these may be unlikely since the point source criteria requires small elongation values. Our sources also do not appear to be close to potential lensing sources (see Figure 2). Moreover, Mason et al. (2015) showed that the effect of magnification bias on the luminosity function determination is small. Shibuya et al. (2022) calculated the merger fraction of 10% to 70% for bright −24 M U V −22 galaxies (Harikane et al. 2022). Considering the inferred M U V of our sources, there is a non-zero possibility of merger contaminants, especially since galaxy formation in the early universe may involve major mergers.
If the point sources are revealed as quasars by future spectroscopic follow-ups, then it may imply substantial population of low-luminous quasars in the early universe. Since quasar activity is associated with rapid accretion, this may allow a pathway for the rapid formation of massive black hole seeds. Distinguishing between compact sources and contaminants may require higher spatial resolution than the capabilities of HST. What is clear is that point sources selected from deep 3D-HST and medium-deep SuperBoRG both consistently suggest a bright end excess. With the derived upper limit in their number density, we also find that the number of the point-source population dominates only ∼ < 5 % over the galaxy population at M UV ∼ > − 21 mag (Figure 7). This suggests that their contribution to the faint end luminosity function, and thus cosmic reionization, is likely limited.

Caveat: low-redshift interlopers
As demonstrated in this paper, the color and morphology selection using current HST capabilities is still challenging in distinguishing between z ∼ 7-8 sources and Galactic interloper stars. The difficulty is further compounded by the fact that these objects are detected at low signal to noise ratio. We also note that the astrometry analysis did not show significant differences in the apparent proper motion following the analysis used in Morishita et al. (2020), so high-quality spectroscopic analysis was the only reliable, albeit incomplete, metric in eliminating low redshift contaminants. We also calculated β U V , defined in Eq.6, of all SpeX dwarf star template spectra, in the event they are misidentified as galaxies. We find a mean slope ofβ U V = −1.2 ± 1.2, which is nearly indistinguishable from those of galaxies . However, the range of predicted β U V is also large, spanning between −5.32 β U V +6.76, which suggests that β U V is not a useful metric to distinguish galaxies (including quasars) from dwarf stars. In fact, the β U V of spectroscopically confirmed dwarf stars (see Appendix C) have a mean slope ofβ U V = −4.1. Thus, we require spectroscopic follow-up to confirm the redshifts and spectroscopic properties of the targets identified in this study. Therefore, for now, our luminosity function estimates serve as an upper limit to the quasar number density at z ∼ 7-8.
Despite this fact, our study shows that low-resolution spectroscopy around 1µm is an effective method in identifying foreground stars. It is noted that when we opt to use the Y 105 /J 125 /H 160 color selection (Morishita et al. 2020;Morishita 2021;Roberts-Borsani et al. 2021) for the field available (i.e. the GOODS-South field), no dwarf star contaminants are confirmed by the grism data (Ishikawa, in prep.). This may be indicative that Y 105 is an effective z > 7 selector in the absence of other filter observations as proposed by Morishita (2021). Despite this challenge, we eliminate a few dwarf stars contaminants with grism spectroscopy. For very faint point sources that do not have sensitive spectroscopic data, we demonstrate that the SED fits favor z ph ∼ 7-8 sources over dwarf stars.

CONCLUSION
We searched for z ∼ 8 Lyman dropout point sources with the archival 3D-HST data. 3D-HST surveys nearly 700 square arcminutes of the CANDELS fields, reaching the depth of J 125 and H 160 ∼ 26 mags. We combined Lyman dropout color and point source selections, with additional photometric redshift estimates and grism spectroscopy to eliminate low-redshift contaminants, and identified three z ∼ 7-8 point source candidates.
We then investigated the physical properties of the point sources by using the available multi-band photometric data. SED analyses suggest that these sources are potentially quasars or compact star-forming galaxies. Assuming these sources are galaxies, the fitting results revealed high star formation surface density. This is consistent with the redshift trend of previously identified luminous galaxies of comparable luminosity, M U V ∼ −22; however, we find even larger Σ SFR values than those predicted by the relation. We measured the EW Hβ+[O iii] emission for the 3 targets with Spitzer /IRAC photometry available and found moderately high equivalent widths of ∼ 500 − 1000Å.
We calculated the number density distribution and derived the luminosity function of the point sources. Similar to previous HST z ∼ 7-8 point source surveys, we found that the inclusion of point sources revealed an excess in the bright end at M U V −22, consistent with the SuperBoRG point source survey. We combined the z ∼ 7-8 point source point source candidates with published galaxy number densities to estimate the total galaxy luminosity function. We found that the bestfitting models all point to a bright M * U V cut-off, which departs from the known galaxy luminosity function.
The deeper observations of 3D-HST allowed us to extend the dynamic range covered by our previous work in SuperBoRG. We did not identify point sources in the faint end. Moreoever, we found that they make up less than 5% of the galaxy fraction at the faint end. Thus it is unlikely that these point sources have major contributions to cosmic reionization.
If these z ∼ 7-8 point sources are confirmed to be quasars, our results suggest that quasars in this luminosity range may be more abundant in the early universe. If they turn out to be a galaxy population, it would indicate the presence of compact and intense star-formation in the early universe. Further follow-up observations are required to confirm the inferred properties of our point sources. Future surveys using the infrared optimized James Webb Space Telescope and large field-of-view capable Roman Space Telescope may resolve the current limitations of HST.
ACKNOWLEDGMENTS Support for this work was provided by NASA through grant numbers HST-GO-15212.002 and HST-AR-15804.002-A from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 5-26555. This work is based on observations taken by the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.
This work is based on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Table 9. The HST and Spitzer /IRAC aperture photometry taken at 0. 7 apertures of the spectroscopically confirmed dwarf stars. These targets meet the color dropout and point source morphological selections. Non-detections for both HST and Spitzer are shown at 2σ upper limits. The uncertainties of real detections are defined at 1σ. The best-fit HST grism spectra, which also correspond to fluxes in F140W, are shown in Figure 10.  Figure 10. The observed 3D-HST grism spectra in black with the best-fit SpeX templates in red.
note that the best fitting templates are approximate and are only used to identify dwarf star contaminants. Accurate characterization will require careful stellar modeling, as demonstrated by Aganze et al. (2022). With the exception of COS 25286, all of these sources were identified by Aganze et al. (2022); it is possible that COS 25286 did not make their selection due to its low S/N ratio. These targets have F125W magnitudes ranging between 25 I 125 22 mags. While three out of four dwarf stars are much brighter than our final targets, some like COS 25286 can appear as faint as z ∼ 8 candidates. If we calculate the β U V of these targets using Eq.6, we find very blue slopes ofβ U V = −4.1 ± 0.5. We show the target fluxes in Table 9; the grism spectra and their best fitting SpeX templates (Rayner et al. 2003;Burgasser 2014) are shown in Figure 10.