Unresolved z ∼ 8 Point Sources and Their Impact on the Bright End of the Galaxy Luminosity Function

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). 3D-HST is an HST near-infrared (NIR) spectroscopic survey designed to study galaxies across the universe. It has surveyed nearly 700 arcmin² 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₁₆₀ reach median 5σ depths at 26 mag at a 1'' aperture. Further details of the survey and the published catalog can be found in Skelton et al. (2014) and Momcheva et al. (2016).

The reason for our choice of using the 3D-HST catalog over the catalogs published by the CANDELS team (Guo et al. 2013; Galametz et al. 2013; Nayyeri et al. 2017; Stefanon et al. 2017b; Barro et al. 2019) is that uniform analysis was performed on all five CANDELS fields by the 3D-HST team to create the source catalog. This vastly simplifies the source detection procedure (Section 2.2) and the completeness simulation analysis (Section 3.3), for calculating the number density of the target population. However, due to inconsistencies in the filter coverage, we only analyze four of the five CANDELS fields (AEGIS, COSMOS, GOODS-South, and UKIDSS-UDS), where the F814W, F125W, and F160W filters are available (Section 2.2). We also exploit the published G141 grism data to identify low-redshift interlopers (Section 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 006 pixel⁻¹. The photometric source catalogs were produced using point-spread function (PSF)-matched aperture photometry, reduced using SExtractor (Bertin & Arnouts 1996), and flux-calibrated to an aperture radius of 07. The HST ACS and WFC3 images were convolved to match the HST/F160W PSF (∼014). Ground-based optical, NIR, and Spitzer/IRAC fluxes were similarly PSF-matched to a combination of F125W, F140W, and F160W priors and aperture-corrected to F160W (or F140W, otherwise).

Our strategy in identifying z ∼ 7–8 point-source candidates from 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) and Morishita (2021) selection, which uses the F105W/F125W/F160W (Y₁₀₅/J₁₂₅/H₁₆₀) filters, the 3D-HST catalog does not include Y₁₀₅ fluxes. Instead, we use the Bouwens et al. (2015) color dropout criteria, which are based on an F814W/F125W/F160W (I₈₁₄/J₁₂₅/H₁₆₀) selection:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{S}}/{{\rm{N}}}_{\mathrm{125,160}}\gt 5.0\\ {\rm{S}}/{{\rm{N}}}_{\mathrm{blue}}\lt 2.0\\ {I}_{814}-{J}_{125}\gt 2.2\\ {J}_{125}-{H}_{160}\lt 0.4\\ {I}_{814}-{J}_{125}\gt 2\times ({J}_{125}-{H}_{160})+2.2.\end{array}\end{eqnarray} \tag{ 1 }$$

Compared to the Morishita et al. (2020) selection, this color selection results in a broader z ∼ 7–8 range. Also, the GOODS-North catalog does not include I₈₁₄ data, so we only examine four of the five CANDELS fields. Where available, we use the HST/ACS blue filters (I₈₁₄ and bluer) to determine the Lyman dropout with strict blue signal-to-noise ratio (S/N) constraints, combined with the nondetection flag from the HST pipeline catalog. Since we require 2σ nondetections for I₈₁₄ fluxes, we calculate the resulting I₈₁₄ − J₁₂₅ lower-limit color as

$$\begin{eqnarray}&&{({I}_{814}-{J}_{125})}_{\mathrm{lim}}=-2.5{\mathrm{log}}_{10}(2{\sigma }_{814}/{J}_{125}).\end{eqnarray} \tag{ 2 }$$

We do not include sources with I₈₁₄ − J₁₂₅ that fall below the nondetection limit (i.e., that fall 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 Section 2.4).

Zoom In Zoom Out Reset image size

From the color dropouts, we identify point sources based on two morphological parameters, elongation and flux concentration, measured in the H₁₆₀ filter. The selection criteria are defined by Morishita et al. (2020):

$$\begin{eqnarray}&&\begin{array}{l}\qquad e\lt 1.2\\ {f}_{5}/{f}_{10}\gt 0.5.\end{array}\end{eqnarray} \tag{ 3 }$$

Elongation, e (the ratio of the semimajor/semiminor axes), describes the circularity of the source, and the flux concentration (the flux ratio between the inner and outer radii) describes the compactness. Morishita (2021) found 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₁₆₀ flux ratios at different radii, we determine the f₅/f₁₀ flux concentration, the flux ratios taken within the 5 pixel (03) and 10 pixel (06) 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 Section 2.4). Although the 3D-HST source catalogs include the star_class flag, which classifies a source as starlike or otherwise, Finkelstein et al. (2015) and Morishita (2021) have demonstrated that this flag is 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., Bouwens et al. 2015; 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₈₁₄/J₁₂₅/H₁₆₀ dropout point sources. Of these point sources, seven meet the photometric redshift selection of z_ph > 7, discussed in Section 2.3. We show the J₁₂₅ − H₁₆₀ versus I₈₁₄ − J₁₂₅ 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 candidate list consists of three point sources, which are listed in Table 1. We also visually inspect the HST images of the z_ph-selected 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₈₁₄ − J₁₂₅ limited, noncandidate point sources in Appendix A.

Zoom In Zoom Out Reset image size

Table 1. Targets Selected from the 3D-HST Fields Based on Color and Morphology Criteria

Target	R.A.	Decl.	e	f5/f10	I814 − J125	J125 − H160	zph	p(zph > 6)	${\chi }_{\nu ,{z}_{\mathrm{ph}}}^{2}$	${\chi }_{\nu ,\mathrm{BD}}^{2}$
	(J2000)	(J2000)	(a/b)	⋯	(AB mag)	(AB mag)	⋯	⋯	⋯	⋯
GDS 45797	53.036877	−27.696156	1.008	0.98 ± 0.01	>3.4	0.31 ± 0.08	${8.1}_{-0.4}^{+0.3}$	99%	4.67	9.42
EGS 515	215.252975	+53.028187	1.111	0.60 ± 0.02	>2.2	−0.40 ± 0.16	${7.5}_{-0.6}^{+0.7}$	99%	2.54	9.70
EGS 29337	215.050369	+53.007496	1.024	0.53 ± 0.02	>2.4	0.02 ± 0.09	${7.5}_{-0.8}^{+0.9}$	99%	1.61	12.45

Note. Photometric redshifts estimated with EAZY. The targets listed here are predicted to have a photometric redshift probability of p(z_ph > 6) > 70%. We also indicate whether grism spectra are 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.

Download table as:  ASCIITypeset image

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 a 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 are excluded in our initial color selection in Section 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 (Section 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) note 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 the 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.

Zoom In Zoom Out Reset image size

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_UV, 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.  Bagpipes SED Fit Priors Assuming a Young Stellar Population Model

Parameter	Units	Priors
z	⋯	EAZY zph
τage	(Gyr)	[0.001, 1]
${\mathrm{log}}_{10}({M}_{\star }/{M}_{\odot })$	⋯	[6, 15]
Z	(Z⊙)	[0, 1]
AV	mag	[0, 1]
${\mathrm{log}}_{10}U$	⋯	−3

Note. The model redshifts are fixed to the EAZY-derived 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, A_V measures the dust attenuation, and U is the nebular ionization parameter, which captures the nebular emission and continuum components.

Download table as:  ASCIITypeset image

Table 4. Best-fit Bagpipes SED Parameters for Each Object

Target	Reff	MUV	βUV	SFR	τage	${\mathrm{log}}_{10}({M}_{* }/{M}_{\odot })$	${\mathrm{log}}_{10}(Z/{Z}_{\odot })$	ΣSFR	EWHβ + [O III ]
	(kpc)	(AB mag)	⋯	(M⊙ yr−1)	(Gyr)	⋯	⋯	(M⊙ yr−1 kpc−2)	(Å)
GDS 45797	<0.7	−22.56	−2.10	12 ± 7	0.5 ± 0.3	9.1 ± 0.2	−1.70 ± 0.43	>34	⋯
EGS 515	0.8 ± 0.7	−21.86	−2.20	15 ± 7	0.6 ± 0.3	9.2 ± 0.3	−0.57 ± 0.43	>9	1200 ± 500
EGS 29337	0.8 ± 0.7	−21.85	−1.40	46 ± 27	0.6 ± 0.2	10.3 ± 0.2	−0.60 ± 0.42	>11	500 ± 200

Note. The model redshifts are fixed to the EAZY z_ph in Table 1. The model priors are listed in Table 3. PSF-limited R_eff are shown as upper limits.

Download table as:  ASCIITypeset image

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 resembles the Lyman break of z ∼ 7–8 objects, making them likely interloper contaminants. There are notable spectral features at 1, 1.25, 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₁₄₀ ≈ 26 mag) or to incomplete spectral coverage (landing on/outside of the detector edge). It is noted that the G141 grism does not cover the redshifted Lyα break at z ∼ 7–8, making 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 Section 2.2, we perform spectral fits to low-mass L and T dwarf template spectra that were observed with the SpeX spectrograph on the NASA Infrared Telescope Facility (Rayner et al. 2003). The template spectra are obtained from the SpeX Prism Library (Burgasser 2014). Based on the spectral fits, we identify four T dwarfs with clear spectral features. Their 1D spectra are shown in Appendix C. In fact, three 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₁₆₀ ≲ 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 the SpeX dwarf star templates using EAZY. In Table 1 we list the χ_ν,BD, and in Figure 3 we show the best-fit SEDs. When these are compared to the χ_ν of the z_ph fits, we find that our final point sources are better constrained as z ∼ 7–8 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.

Postage-stamp images of our final point sources are shown in Figure 2. The images are extracted from the 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) data that make up the 3D-HST catalogs. From the images, we clearly see 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₈₁₄ fluxes, despite meeting the nondetection criteria. After comparing the different filter images, we conclude 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 suggests spurious ground-based blue fluxes observed with Subaru (Taniguchi et al. 2007) for GDS 45797. However, careful examination of the Subaru images suggests that they are artifacts due to diffraction spikes from a nearby star. So, we treat these blue fluxes as nondetections in our analysis. Fortunately, the inclusion of these blue fluxes does not have a major effect on the z_ph or the SED fit results. Another uncertainty comes from the 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 07 apertures. As a result, we individually inspect the IRAC images to decide whether to include them in our analysis. If we visually confirm obvious contamination within a channel, we omit its flux.

Finally, we also cross-match our HST-selected sources 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 targeted observations, at shallower depths, have detected z ∼ 7 low-luminosity quasars (e.g., Bañados et al. 2018a; Wang et al. 2021a); however, no clear X-ray emission was detected at 0.5–7 keV. We place upper limits on the fluxes and luminosities of our targets in Table 5. Moreover, there are known Compton-thick z > 7 quasar candidates with faint X-rays (Fujimoto et al. 2022), so this is not entirely unusual. Detailed X-ray analyses will be left for future study.

We cross-match our point-source selection with findings from Bouwens et al. (2015), which examined the CANDELS, HUDF09, HUDF12, ERS, and BoRG/HIPPIES fields to estimate the galaxy UV luminosity function. Both studies apply similar Lyman-dropout color selections. The main difference is that we explicitly search for point sources with the morphology selection defined in Section 2.2 using the f₅/f₁₀ flux ratios, while Bouwens et al. (2015) implemented the stellarity parameter (i.e., the 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 a separation of δ r = 011, which is well within the PSF uncertainty. In fact, this object has also been spectroscopically confirmed as a galaxy at z_ph = 7.477 (Roberts-Borsani et al. 2016; Stark et al. 2017). This proves 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 conclusions. We note that the other sources were not identified by other z ∼ 7–8 studies (e.g., Bouwens et al. 2015; Roberts-Borsani et al. 2016; Stark et al. 2017). In fact, these sources have larger f₅/f₁₀ values compared to EGS 29337, which suggest that they appear more pointlike, and thus are 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 were often excluded in earlier studies of high-redshift galaxies.

We list the best-fit stellar population properties of the point sources, which are estimated with Bagpipes SED fits, in Table 4. Due to the limitation of the data, we assume that our detected point sources are well represented by a young stellar spectrum. The fit results predict subsolar 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 on the star formation history and age–metallicity–dust degeneracies (e.g., Figure 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, 014, we find that two 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{SFR}}=\displaystyle \frac{\mathrm{SFR}}{2\pi {R}_{\mathrm{eff}}^{2}}.\end{eqnarray} \tag{ 4 }$$

The calculated Σ_SFR serves as a lower limit since its uncertainty is dependent on the upper-limit uncertainty in R_eff. We compare it with the inferred Σ_SFR, which is calculated from the M_UV–SFR relation (Kennicutt 1998; Ono et al. 2013) defined as follows:

$$\begin{eqnarray}&&{M}_{\mathrm{UV}}=-2.5{\mathrm{log}}_{10}\left[\displaystyle \frac{{{\rm{\Sigma }}}_{\mathrm{SFR}}\cdot \pi {R}_{\mathrm{eff}}^{2}}{2.8\times {10}^{-28}({M}_{\odot }\,{\mathrm{yr}}^{-1})}\right]+51.59.\end{eqnarray} \tag{ 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_UV. We also compare the Σ_SFR redshift evolution. Our results appear to be consistent with the trends of high-redshift galaxies discussed in Ono et al. (2013) and Holwerda et al. (2018), which predict a greater number of compact galaxies with high SFRs at earlier epochs.

Using the best-fit Bagpipes SEDs, we calculate the UV continuum slope β_UV. We adopt the formula defined as follows (e.g., Dunlop et al. 2013):

$$\begin{eqnarray}&&{\beta }_{\mathrm{UV}}=-2.0+4.39\times ({J}_{125}-{H}_{160}),\end{eqnarray} \tag{ 6 }$$

where J₁₂₅ and H₁₆₀ are the best-fit magnitudes from the best-fit Bagpipes SEDs. The calculated β_UV are listed in Table 4 with a mean slope of $\bar{{\beta }_{\mathrm{UV}}}=-1.90\pm 0.35$. The resulting $\bar{{\beta }_{\mathrm{UV}}}$ is consistent with the β_UV 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 (EW) due to the Hβ + [O III] emission lines from the best-fit 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 EWs as

$$\begin{eqnarray}&&{\mathrm{EW}}_{{\rm{H}}\beta +[{\rm{O}}\,\mathrm{III}]}\,=\frac{({f}_{\mathrm{ch}2}-{f}_{\mathrm{cont}})}{{f}_{\mathrm{cont}}}\frac{{\rm{\Delta }}{\lambda }_{\mathrm{ch}2}}{(1+{z}_{\mathrm{ph}})},\end{eqnarray} \tag{ 7 }$$

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 FWHM of the CH2 filter, and z_ph is from EAZY. Comparing the values in Table 4 and the SED plots in Figure 3, we see that the EW_{Hβ + [O III ]} estimates can be applied to three targets. These objects show a moderate EW_{Hβ + [O III ]} of 500–1000 Å, similar to estimates from other high-redshift surveys (Labbé et al. 2013; Schenker et al. 2013; Smit et al. 2014, 2015; Morishita et al. 2020). However, this is based on the assumption that z_ph is correct. We compare the redshift evolution of the measured EW_{Hβ + [O III ]} in Figure 4, and our results appear to be consistent with other survey results. Future infrared observations with higher resolution and sensitivity are needed to better characterize the predicted Hβ + [O III] emission.

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 3D-HST. We follow the completeness simulation treatment in Leethochawalit et al. (2021) to calculate V_eff (Oesch et al. 2012; Calvi et al. 2016; Carrasco et al. 2018; Morishita et al. 2018).

We inject 500 sources into each 3D-HST image at each (M_UV, z) bin: 100 ΔM_UV bins across −26 ≤ M_UV ≤ −16 and 13 Δz bins across 7 ≤ z ≤ 9.4. All simulated sources in a given (M_UV, z) grid have the same UV slope, which are randomly drawn from a Gaussian distribution with a mean slope of $\bar{{\beta }_{\mathrm{UV}}}=-2.2\pm 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 functions of the extraction completeness in Figure 5.

Zoom In Zoom Out Reset image size

After we calculate V_eff, we estimate the number density of the point sources shown in Table 6 within each ΔM_UV = 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 it with results from previous surveys of point sources and galaxies at z ∼ 7–8.

Zoom In Zoom Out Reset image size

We fit the point-source number density with both the Schechter function (Equation (8); Schechter 1976) and a double power law (DPL) (Equation (9); Hopkins et al. 2007), where ϕ* is the characteristic normalization, ${M}_{\mathrm{UV}}^{* }$ is the characteristic UV luminosity defined at M_UV (1450 Å), α defines the faint-end slope, and β defines the bright-end slope:

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{\mathrm{Sch}} & = & \displaystyle \frac{\mathrm{ln}10}{2.5}\phi * \times {10}^{-0.4(M^{\prime} -{M}_{* })(\alpha +1)}\\ & & \times \,\exp \left[-{10}^{-0.4({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}\right]\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{\mathrm{DPL}} & = & \displaystyle \frac{\mathrm{ln}10}{2.5}{\phi }^{* }\times \left[{10}^{0.4(\alpha +1)({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}\right.\\ & & \left.+\,{10}^{0.4(\beta +1)({M}_{\mathrm{UV}}-{M}_{\mathrm{UV}}^{* })}\right].\end{array}\end{eqnarray} \tag{ 9 }$$

To properly include the bins of nondetection of the point-source population at the faintest and brightest ends in fitting evaluation, we incorporate the upper limits from nondetections following the derivations from Sawicki (2012) for χ² minimization. The derivations for M_UV and the χ² minimization are shown in Appendix D. We perform Markov Chain Monte Carlo sampling, using the 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}_{\mathrm{UV}}^{* }$. 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 DPL model fits, we fix the faint-end slope to α = −1.2 and the bright-end slope to β = −2.7 based on extrapolations of 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 luminosity function parameters of the point sources are shown in Table 7. If we freely fit for M_UV, we find that the DPL fit produces more reasonable parameters (${M}_{\mathrm{UV}}^{* }\approx -24$) than the best-fit Schechter function, which instead suggests an unrealistically bright UV cutoff (${M}_{\mathrm{UV}}^{* }\approx -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_UV fit, then it is difficult to confidently favor either function 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).

Table 7. Best-fit Luminosity Function Parameters

Survey	Model	ϕ*	${M}_{\mathrm{UV}}^{* }$	α	β
		(Mpc−3 mag−1)	(AB mag)
Point sources: this work	Schechter	$({1.5}_{-1.0}^{+4.3})\times {10}^{-8}$	$-{24.8}_{-0.9}^{+1.1}$	[−2.2]	⋯
	DPL	$({3.7}_{-1.7}^{+4.3})\times {10}^{-7}$	$-{23.5}_{-1.4}^{+1.2}$	[−1.2]	[−2.7]
Point sources: this work + SuperBoRG	Schechter	$({1.0}_{-0.7}^{+3.4})\times {10}^{-8}$	$-{24.7}_{-0.9}^{+1.2}$	[−2.2]	⋯
	DPL	$({2.1}_{-1.1}^{+3.3})\times {10}^{-7}$	$-{23.6}_{-1.5}^{+1.5}$	[−1.2]	[−2.7]
Point sources + galaxies (fixed α, β)	Schechter	$({2.7}_{-0.8}^{+1.0})\times {10}^{-5}$	$-{21.9}_{-0.2}^{+0.2}$	[−2.2]	⋯
	DPL	$({7.7}_{-0.3}^{+0.4})\times {10}^{-3}$	$-{17.6}_{-0.2}^{+0.2}$	[−1.2]	[−2.7]
Point sources + galaxies (freely fit)	Schechter	$({3.4}_{-2.9}^{+6.3})\times {10}^{-6}$	$-{22.8}_{-1.0}^{+0.5}$	$-{2.5}_{-0.1}^{+0.1}$	⋯
	DPL	$({4.1}_{-3.1}^{+3.8})\times {10}^{-4}$	$-{20.0}_{-0.8}^{+0.4}$	$-{2.0}_{-0.3}^{+0.3}$	$-{3.6}_{-0.3}^{+0.3}$

Note. For the point sources, we fix both the α and β slopes depending on the model used and freely fit for ϕ* and ${M}_{\mathrm{UV}}^{* }$. Parameters that are fixed are shown in square brackets. The Schechter model slope is fixed to the galaxy faint-end slope from Bouwens et al. (2021), and the DPL slopes are fixed to the active galactic nucleus 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.

Download table as:  ASCIITypeset image

Finally, we fit the combined point-source and galaxy luminosity functions. Here, we remove EGS 29337 from the point-source luminosity function since it is already included in the Bouwens et al. (2015) determination (discussed in Section 3.1). First, we fit with the slopes fixed, and then we freely fit over all parameters. With the current survey volume by HST, it is difficult to confidently favor either function over the other for the combined luminosity function. With the Schechter model, both the fixed-slope and freely fit runs produce brighter UV cutoffs at ${M}_{\mathrm{UV}}^{* }\approx -21.9$, compared to ${M}_{\mathrm{UV}}^{* }=-22.8$ by Bouwens et al. (2021), suggesting an excess of bright sources. The freely fit DPL model produces a superposition of point-source Schechter and galaxy Schechter functions also with a slight excess of ${M}_{\mathrm{UV}}^{* }=-20.0$ and a steep β = −3.6.

The exact shape of the galaxy luminosity function is under debate. For example, Harikane et al. (2022) performed a two-component luminosity function fit (i.e., DPL+DPL or DPL+Schechter) to the combined quasar+galaxy populations. In contrast, we model a single-component function across all magnitudes for both the point-source and galaxy number densities (i.e., Schechter or DPL). The main difference between the two studies is that Harikane et al. (2022) extrapolated 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_UV ≲ −24 SuperBoRG 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare the point-source number density and luminosity function against the galaxy luminosity function measured at z ∼ 7–8 (Bouwens et al. 2021). We quantify the fraction of extended galaxies among all sources (galaxies and point sources) as a function of M_UV as follows:

$$\begin{eqnarray}&&{f}_{\mathrm{galaxy}}({M}_{\mathrm{UV}})=\displaystyle \frac{{\phi }_{\mathrm{galaxy}}}{{\phi }_{\mathrm{galaxy}}+{\phi }_{\mathrm{points}}}\end{eqnarray} \tag{ 10 }$$

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 them with the results from Harikane et al. (2022). We find that these point sources dominate at the bright M_UV magnitudes. This suggests that the bright-end excess implied 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.

With the exception of EGS 29337, it is currently difficult to distinguish whether our final candidates are nonactive galaxies or quasar-hosting galaxies, as the inferred sizes of the point sources are 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 nonactive galaxies, we predict a high Σ_SFR based on our SED fitting analysis. Previously, Oesch et al. (2010) and Ono et al. (2013) observed constant Σ_SFR from z ∼ 4 to z ∼ 7 with a weak increase toward higher redshifts. Our SED analysis of the point sources appears to support this increasing Σ_SFR trend, although we predict even larger values, as shown in Figure 4. Finally, given their predicted SFRs and stellar masses, these sources may be progenitors of massive quiescent galaxies that were 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 galaxies with a high SFR (Roberts-Borsani et al. 2016; Stark et al. 2017).

While much of our SED analysis assumes the star-forming 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, 2019; Inayoshi et al. 2020). Either there is an enhancement in the population of unobscured quasars or some physical mechanisms, such as powerful outflows, may drive obscured quasars to appear as luminous as unobscured quasars. Also, while nondetections 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 the z ∼ 7 sources identified by Laporte et al. (2017). In fact, the possibility of either a compact starburst or a quasar is supported by the recent discovery of a 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 to reach below the predicted sizes and magnitudes of our sources at the range of −22 ≲ M_UV ≲ −18. Marshall et al. (2021) predicted that deep observations by NIRCam may allow us to study the quasars and their host galaxies at z ∼ 7. Its imaging and spectroscopic capability may enable us to confidently distinguish the sources as quasars or as compact star-forming galaxies. If they are revealed as quasars, they will become some of the most distant quasars ever discovered.

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 galaxy surveys. This may support the existence of a bright-end excess in the early universe. While ground-based observations (e.g., Bowler et al. 2020) 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) presented a comprehensive analysis of the luminosity function by combining the quasar and galaxy populations identified in the Hyper Suprime-Cam 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 activity 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 of compact star-forming galaxies or quasars are consistent with these possibilities.

Other possibilities include the superposition of lensed galaxies or even merging galaxies. However, these may be unlikely since the point-source criteria require 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 luminosity function determination is small. Shibuya et al. (2022) calculated a merger fraction of 10% to 70% for bright −24 ≲ M_UV ≲ −22 galaxies (Harikane et al. 2022). Considering the inferred M_UV of our sources, there is a nonzero possibility of merger contaminants, especially since galaxy formation in the early universe may involve major mergers.

If the point sources are revealed to be quasars by future spectroscopic follow-ups, then a substantial population of low-luminosity quasars may be inferred to exist 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 both deep 3D-HST and medium–deep SuperBoRG consistently suggest a bright-end excess. With the derived upper limit in their number density, we also find that the number density 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 to cosmic reionization, is likely limited.

As demonstrated in this paper, color and morphology selection using current HST capabilities is still challenging when 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 S/N. We also note that the astrometry analysis does not show significant differences in the apparent proper motion following the analysis used in Morishita et al. (2020), so high-quality spectroscopic analysis is the only reliable, albeit incomplete, metric in eliminating low-redshift contaminants. We also calculate the β_UV, defined in Equation (6), of all SpeX dwarf star template spectra, in case the stars are misidentified as galaxies. We find a mean slope of $\bar{{\beta }_{\mathrm{UV}}}=-1.2\pm 1.2$, which is nearly indistinguishable from those of galaxies (Dunlop et al. 2013). However, the range of predicted β_UV is also large, spanning −5.32 ≲ β_UV ≲ +6.76, which suggests that β_UV is not a useful metric to distinguish galaxies (including quasars) from dwarf stars. In fact, the β_UV of spectroscopically confirmed dwarf stars (see Appendix C) have a mean slope of $\bar{{\beta }_{\mathrm{UV}}}=-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 for identifying foreground stars. It is noted that when we opt to use the Y₁₀₅/J₁₂₅/H₁₆₀ 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 (Y. Ishikawa 2022, in preparation). This may be indicative of the observation that Y₁₀₅ 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 star 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.