A Search for H-Dropout Lyman Break Galaxies at z ∼ 12–16

Observing the first galaxy formation is one of the main goals in modern astronomy. One of the most straightforward approaches to achieve this goal is to observe forming galaxies directly in the early universe. Large telescopes currently in operation have yielded the most distant objects so far. These highest-redshift objects have posed various interesting questions for astronomy. For example, the most distant quasars at z > 7 raised a question about how to form black holes as massive as ∼10⁹ M_⊙ in the limited cosmic time (e.g., Mortlock et al. 2011; Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021). Thus, searching for the most distant objects not only is the frontier of the knowledge of human beings but also has great power to reveal the formation physics of various objects in the early universe (e.g., see review by Stark 2016; Dayal & Ferrara 2018; Robertson 2021).

The current record of the highest-redshift galaxy spectroscopically confirmed is GN-z11 at z ∼ 11 measured with detections of the Lyman break and rest-frame ultraviolet (UV) metal lines (Oesch et al. 2016; Jiang et al. 2021). A major surprise of GN-z11 is its remarkably high luminosity, M_UV = −22.1 mag. Given that it is not gravitationally lensed, GN-z11 is located in the brightest part of the rest-frame UV luminosity function. Although the narrow field of view (FOV) of Hubble Space Telescope (HST)/Wide Field Camera 3 (WFC3) in the near-infrared has limited the imaging survey areas to <1 deg², several studies using HST report very luminous Lyman break galaxies (LBGs) at z ∼ 9–10 more frequently than the expectation from a Schechter-shape luminosity function (e.g., Morishita et al. 2018; Finkelstein et al. 2021a; see also Roberts-Borsani et al. 2021a). More statistically robust results have come from few-square-degree near-infrared imaging surveys with the Visible and Infrared Survey Telescope for Astronomy (VISTA) and UK Infrared Telescope (UKIRT) such as UltraVISTA (McCracken et al. 2012), the UKIRT InfraRed Deep Sky Surveys (UKIDSS; Lawrence et al. 2007), and the VISTA Deep Extragalactic Observation (VIDEO) survey (Jarvis et al. 2013). These surveys have revealed that the UV luminosity functions at z ∼ 9–10 are more consistent with a double-power-law function than a standard Schechter function (Stefanon et al. 2017, 2019; Bowler et al. 2020). Previous studies also report similar number density excesses beyond the Schechter function at z ∼ 4–7 (Ono et al. 2018; Stevans et al. 2018; Adams et al. 2020; Harikane et al. 2022), implying little evolution of the number density of bright galaxies at z ∼ 4–10 (Bowler et al. 2020; Harikane et al. 2022). Although spectroscopic observations are required to confirm these results, the studies indicate that there are a larger number of luminous galaxies at z ∼ 9–11 than previously thought, which formed in the early universe of z > 10.

In addition to these observations of bright galaxies at z ∼ 9–11, several studies independently suggest the presence of star-forming galaxies in the early universe even at z ∼ 15. A candidate for a z ∼ 12 galaxy is photometrically identified in very deep HST/WFC3 images obtained in the Hubble Ultra Deep Field 2012 (UDF12) campaign (Ellis et al. 2013). Balmer breaks identified in z = 9–10 galaxies indicate mature stellar populations whose ages are ∼300–500 Myr, implying early star formation at z ∼ 14–15 (Hashimoto et al. 2018; Laporte et al. 2021; see also Roberts-Borsani et al. 2020). An analysis of passive galaxy candidates at z ∼ 6 reports that their stellar population is dominated by old stars with ages of ≳700 Myr, consistent with star formation activity at z > 14 (Mawatari et al. 2020a).

Motivated by these recent works, we search for H-band dropout (H-dropout) LBGs whose plausible redshifts are z ∼ 12–16. ¹⁶ Given the observed number density of luminous galaxies at z ∼ 9–10 and its little redshift evolution from z ∼ 4 to 10, it is possible that one to several z ∼ 12–16 galaxies will be found in currently available data sets obtained by surveys with large ground- and space-based telescopes. This search for z ∼ 12–16 galaxies is important not only for understanding early galaxy formation but also for designing survey strategies with upcoming space missions that will study the z > 10 universe such as the James Webb Space Telescope (JWST).

This paper is organized as follows. We describe photometric data sets and a selection of z ∼ 12–16 galaxies in Section 2 and ALMA follow-up observations for one of our candidates in Section 3. Results of spectral energy distribution (SED) fitting and the UV luminosity function are presented in Sections 4 and 5, respectively. We discuss future prospects with space missions based on our results in Section 6, and we summarize our findings in Section 7. Throughout this paper, we use the Planck cosmological parameter sets of the TT, TE, and EE+lowP+lensing+ext results (Planck Collaboration et al. 2016): Ω_m = 0.3089, Ω_Λ = 0.6911, Ω_b = 0.049, h = 0.6774, and σ₈ = 0.8159. All magnitudes are in the AB system (Oke & Gunn 1983).

2.1. Data Set

2.2. Selection of H-Dropout Galaxies

We construct multiband photometric catalogs in the COSMOS and SXDS fields. We start from K _s (K)-band detection catalogs made by the UltraVISTA (UKIDSS) team using SExtractor (Bertin & Arnouts 1996). We select sources detected in the K _s (K) bands at >5σ levels and not detected in the J band at >2σ levels in a 2''-diameter circular aperture. Then, we measure magnitudes of these sources in the grizyJHK _s (K) images using a 2''-diameter circular aperture centered at their coordinates in the catalogs. Since point-spread functions (PSFs) of the Spitzer/IRAC [3.6] and [4.5] images are relatively large (∼17), source confusion and blending are significant for some sources. To remove the effects of the neighbor sources on the photometry, we first generate residual IRAC images where only the sources under analysis are left by using T-PHOT (Merlin et al. 2016), in the same manner as Harikane et al. (2018, 2019). As high-resolution prior images in the T-PHOT run, we use HSC grizy stacked images whose PSF is ∼07. Then, we measure magnitudes in the IRAC images by using a 3''-diameter apertures in the same manner as Harikane et al. (2018). To account for the flux falling outside the aperture, we apply aperture corrections derived from samples of isolated point souces in Harikane et al. (2018).

To search for z ∼ 12–16 galaxies whose Lyman breaks are redshifted to ∼1.6–2.1 μm, we select H-dropout LBGs from the multiband photometric catalogs constructed above. We adopt the following color selection criteria in the COSMOS and SXDS fields, respectively:

$$\begin{eqnarray}&&H-{K}_{{\rm{s}}}\gt 1.0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{K}_{{\rm{s}}}-[3.6]\lt 0.1,\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&H-K\gt 1.0,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&K-[3.6]\lt 0.1.\end{eqnarray} \tag{ 4 }$$

As shown in Figure 1, these color criteria can select sources at z > 12 while avoiding color tracks of z = 0 − 7 galaxies and stellar sources. To remove foreground interlopers, we exclude sources with detections at >2σ levels in the grizyJ band images. Note that we use the same values for the color criteria in the COSMOS and SXDS fields. Since the filter response profiles are different in the VISTA JHK _s filters for the COSMOS field and the UKIRT JHK filters for the SXDS field, the selection functions will not be identical. We will account for this difference by separately evaluating the selection functions in the COSMOS and SXDS fields in Section 5.1.

Zoom In Zoom Out Reset image size

To remove foreground interlopers further, we conduct a photometric redshift analysis using BEAGLE (Chevallard & Charlot 2016). We adopt a constant star formation history (SFH) with the Chabrier (2003) initial mass function (IMF); stellar ages of 10⁶, 10⁷, and 10⁸ yr; metallicities of 0.2 and 1 Z_⊙; and the Calzetti et al. (2000) dust attenuation law with the V-band optical depth of τ _V = 0−2 (steps of 0.2) to account for very dusty low-redshift interlopers. We select objects whose high-redshift solution is more likely than the low-redshift ones at a >2σ level, corresponding to Δχ² > 4.0, in the same manner as Bowler et al. (2020). Then, we visually inspect images and SEDs of the selected sources to remove spurious sources, sources affected by bad residual features in the T-PHOT-made IRAC images, and extremely red sources (e.g., K _s (K)−[4.5] ≳ 1) that are not likely to be z ∼ 12–16 galaxies.

After these careful screening processes, we finally identify two z ∼ 12–16 galaxy candidates, HD1 and HD2, in the COSMOS and SXDS fields, respectively. Figures 2 and 3 show images and SEDs of HD1 and HD2, and Table 2 summarizes their measured fluxes. HD1 and HD2 are spatially isolated from other nearby sources, ensuring the robustness of the photometry.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Photometry of Our H-Dropout Galaxy Candidates

Name	R.A.	Decl.	Subaru						VISTA/UKIRT				Spitzer
			g	r	i	z	y		J	H	K s /K		[3.6]	[4.5]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)		(9)	(10)	(11)		(12)	(13)
HD1	10:01:51.31	02:32:50.0	<27	<32	<46	<63	<157		<107	<145	510 ± 93		531 ± 108	494 ± 136
HD2	02:18:52.44	−05:08:36.1	<25	<33	<46	<68	<133		<95	296 ± 76	821 ± 63		888 ± 88	1252 ± 132

Note. Col. (1): name. Col. (2): right ascension. Col. (3): decl. Cols. (4)–(13): flux densities in nJy or 2σ upper limits.

Download table as:  ASCIITypeset image

HD1 is also found in the COSMOS2020 catalog (Weaver et al. 2022). However, the photometric redshift of HD1 is 3.6 in the COSMOS2020 catalog. This is due to the difference in the measured magnitudes in the IRAC images. Our measured magnitudes are 24.6 and 24.7 mag in the [3.6] and [4.5] images, respectively, while 24.2 and 23.9 mag (24.4 and 24.1 mag) are cataloged in the COSMOS2020 CLASSIC (FARMER) catalog with very small flux errors of 3−12 nJy. If we remeasure magnitudes by using a larger aperture in the original IRAC images before the T-PHOT run, the magnitudes become brighter owing to a neighboring source located ∼35 from HD1. We additionally test with deeper SMUVS images (Ashby et al. 2018). We carefully measure the magnitudes in the SMUVS images and still find that the best photometric redshift for HD1 is z > 12. Magnitudes in the other bands, including the K _s band, in the COSMOS2020 catalog are consistent with our measurements, although their flux errors are much smaller than ours. Thus, in this study we adopt our measured magnitudes for HD1.

HD1 and HD2 will be observed in a JWST program (GO-1740; Harikane et al. 2021). In addition to these two candidates, the program will target another source, HD3 (R.A. = 02:16:54.48, decl. = −05:09:37.1), which is also a good candidate for a z ∼ 12–16 galaxy with its prominent break with H−K > 1.2 and the best photometric redshift of z_phot = 14.6. However, due to its relatively red color with K−[3.6] = 0.2 ± 0.2, HD3 is not included in our final sample in this paper.

We observed one of the candidates, HD1, in an ALMA Director's Discretionary Time (DDT) program (2019.A.00015.S; PI: A. K. Inoue). Following successful detections of O iii 88 μm emission lines in high-redshift galaxies (e.g., Inoue et al. 2016; Carniani et al. 2017; Laporte et al. 2017, 2021; Marrone et al. 2018; Hashimoto et al. 2018, 2019; Walter et al. 2018; Tamura et al. 2019; Harikane et al. 2020; see also Inoue et al. 2014b), we conducted a spectral scan targeting the O iii 88 μm line using four tuning setups with Band 6 covering the redshift range of 12.6 < z < 14.3. The antenna configurations were C43-2, C43-3, C43-4, and C43-5, and the typical beam size is ∼04–07. We used four spectral windows with 1.875 GHz bandwidths in the Frequency Division Mode and the total bandwidth of 7.5 GHz in one tuning setup. The velocity resolution was set to ∼10 km s⁻¹. The data were reduced and calibrated using the Common Astronomy Software (CASA; McMullin et al. 2007) pipeline version 5.4.0 in the general manner with scripts provided by the ALMA observatory.

Figure 4 shows the obtained spectrum for HD1 extracted from a 10-radius circular aperture. Although there is no signal at a >5σ level, we find a 4σ tentative line-like feature around 238 GHz. As shown in the top panel of Figure 5, this feature is at 237.8 GHz, and the significance level of the peak intensity is 3.8σ in the moment 0 map shown in the bottom panel. Although there are some other line-like features (e.g., 246.3 GHz), the feature at 237.8 GHz has the highest signal-to-noise ratio among the ones in the frequencies free from severe atmospheric O₃ absorption. If this feature is the O iii 88 μm emission line, the redshift of HD1 is z = 13.27, in good agreement with the photometric redshift estimate. A relatively broad line width (∼400 km s⁻¹ in FWHM) is in fact comparable to similarly bright LBGs at z ∼ 6 (Harikane et al. 2020). The emission feature is cospatial with the rest-frame UV emission in the K _s -band image (bottom panel of Figure 5). The integrated line flux is 0.24 ± 0.06 Jy km s⁻¹ or (1.9 ± 0.5) × 10⁻¹⁸ erg s⁻¹ cm⁻², and the line luminosity is L_[OIII] ≃ 3.3 × 10⁸ L_⊙ if z = 13.27 is assumed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The line luminosity is very small compared to the UV luminosity. Since the UV luminosity of HD1 is L_UV =4.8 × 10¹¹ L_⊙, the [O iii]-to-UV luminosity ratio is L_[OIII]/L_UV ∼ 7 × 10⁻⁴. This ratio is the smallest among the galaxies observed in [O iii] 88 μm emissions in the reionization epoch, as well as in the local universe so far (e.g., Inoue et al. 2016; Binggeli et al. 2021). Since the [O iii]-to-UV ratio depends on the oxygen abundance (Harikane et al. 2020), this low ratio indicates a metallicity as low as ∼0.01–0.1 Z_⊙.

Another possibility is that the ALMA line scan just missed the true emission line and the redshift is out of the range of 12.6 < z < 14.3. As we will see in Section 5.1, the redshift selection function is as broad as 12 < z < 17. For HD1, the lower-redshift case (z < 12.6) is not very favored by the SED fitting, but the higher-redshift case (z > 14.3) is still equally likely, as we discuss later in Section 4. Therefore, additional spectroscopic data are highly desired to confirm redshifts of HD1 and HD2. We plan to conduct follow-up observations for the tentative signal in HD1 and to newly obtain spectroscopic data for HD2 in ALMA cycle 8 (2021.1.00207.S; PI: Y. Harikane). We will also observe these candidates with JWST (GO-1740; Harikane et al. 2021), which allows us to examine a wider redshift range than ALMA.

The dust continuum of HD1 remains nondetection, which is consistent with the low-metallicity interpretation from the low L_[OIII]/L_UV ratio discussed above. The obtained 1σ noise level is 8 μJy beam⁻¹. Assuming that HD1 is not resolved in this observation, we obtain a 3σ upper limit on the dust continuum of <24 μJy.

To examine the photometric redshifts of HD1 and HD2 more carefully, we perform a comprehensive SED fitting analysis from optical to submillimeter wavelength using PANHIT (Mawatari et al. 2020b). PANHIT takes the energy conservation of the dust absorption in the rest-frame UV to near-infrared range and the emission in the far-infrared to submillimeter range into account. PANHIT deals with the upper limits for nondetection bands, following the probability distribution function formula proposed by Sawicki (2012). We adopt 1σ for the upper bound of the integral of the probability distribution. In addition to the fluxes in grizyJHK _s (K)[3.6][4.5] measured in Section 2.2, we utilize far-infrared and submillimeter data of the Herschel survey (Oliver et al. 2012) of HD1 and HD2 and the ALMA data obtained for HD1 (see Section 3). Since dust continua of HD1 and HD2 are not detected in these data, we use the upper limits for the SED fitting.

We assume a delayed-τ model for the SFH covering a wide range of histories, including a short-timescale burst, rising, declining, and almost constant cases (Speagle et al. 2014). It is important to include passive galaxy models because the red H−K _s (K) color can be produced by the Balmer break as well as the Lyman break. This may be a major contamination case in our H-dropout selection. Template spectra include the BC03 stellar population synthesis model (Bruzual & Charlot 2003) with the Chabrier (2003) IMF of 0.1–100 M_⊙, the nebular continuum and line emission model (Inoue 2011), and the dust thermal emission with a modified blackbody function. The dust temperature is assumed to be 30 K, 50 K, or 80 K to account for possibilities of dusty interlopers with various temperatures, and the dust emissivity index is fixed at β_dust = 2.0. The effect of the cosmic microwave background on the dust emission (da Cunha et al. 2013) is also taken into account.

The considered fitting parameters are as follows: the SFH timescale is τ_SFH = 0.01, 0.03, 0.06, 0.1, 0.3, 0.6, 1, 3, 6, and 10 Gyr (10 cases); the metallicity is Z = 0.0001, 0.0004, 0.004, 0.008, 0.02 (=Z_⊙), and 0.05 (six cases); the dust attenuation is A _V = 0.01–10 with 20 logarithmic steps; the stellar population age is 7 cases in 1–10 Myr, 8 cases in 10–100 Myr, 15 cases in 100 Myr–1 Gyr, and 8 cases in 1–15 Gyr, but limited by the cosmic age at the redshift of interest; and the redshift is 0.1–20.0 with a 0.1 step assuming a flat prior.

Figure 2 shows the results of the SED fitting analyses, and Table 3 summarizes the results. The best photometric redshifts are always z > 12 for both HD1 and HD2 thanks to the sharp discontinuity between H and K _s (K) bands. The low-redshift solutions are found at z ∼ 4 for both objects with larger χ² values than the z > 12 solutions. These are Balmer break galaxy solutions, and the dust temperature does not affect them because these solutions have very weak or no dust emission. Another type of possible solution are dusty Hα emitters at z ∼ 2, although these solutions are not supported by the nondetections in the far-infrared and submillimeter bands. In these solutions, a strong Hα line boosts K _s (K) band and makes H−K _s (K) color as red as z ∼ 12–16 galaxies. In the very high dust temperature case of 80 K, this solution gives a slightly smaller χ² than those of the Balmer break ones at z ∼ 4, while the lower, more normal dust temperature cases do not favor this type of solution. Moreover, even the 80 K case is significantly less likely compared to the solutions at z > 12 (Δχ² > 4).

Table 3. Physical Properties of Our H-Dropout Candidates

Name	zphot	MUV	SFRUV	$\mathrm{log}{M}_{* }$	A V
(1)	(2)	(3)	(4)	(5)	(6)
HD1	${15.2}_{-2.1(2.7)}^{+1.2(1.6)}$ a	−23.3 a	110 a	∼9–11	<0.08
HD2	${12.3}_{-0.3(0.7)}^{+0.4(0.7)}$	−23.8	170	∼9.8–11	≲0.8

Notes. Col. (1): name. Col. (2): the best photometric redshift with 1σ (2σ) errors. Col. (3): absolute UV magnitude in units of mag. Col. (4): SFR estimated from the UV magnitude by using Equation (7) in units of M_⊙ yr⁻¹. Cols. (5) and (6): stellar mass and dust attenuation suggested by the SED fitting in units of M_⊙ and mag, respectively. See Section 4 for details.

^a z = 13.27 is suggested by the ALMA observations for HD1 (see Section 3), consistent with the photometric redshift estimate within 1σ. The absolute UV magnitude and SFR in this table are calculated based on the assumption of z = 13.27.

Download table as:  ASCIITypeset image

For HD1, the best-fit redshift is z = 15.2 (χ² = 4.7), which is in fact out of the ALMA [O iii] 88 μm scan (Section 3). The case of z = 13.3, corresponding to the possible line feature at z = 13.27, gives χ² = 5.4. Since it is roughly equally likely (within the 1σ confidence range; see Table 3), we show this case in Figure 2. The physical properties are not well constrained, except for the dust attenuation that is A _V < 0.08 (2σ). The stellar mass (M_*) is (1–100) × 10⁹ M_⊙, depending on the stellar age that is not constrained. When the age is less than ∼10 Myr, the stellar mass and star formation rate (SFR) are estimated to be M_* ∼ 1 × 10⁹ M_⊙ and SFR ∼ 10²−10³ M_⊙ yr⁻¹, respectively. For an age of 10–100 Myr (>100 Myr), M_* ∼ (1−10) × 10⁹ M_⊙ and SFR ∼ 10² M_⊙ yr⁻¹ (M_* ∼ (10−100) × 10⁹ M_⊙ and SFR < 10² M_⊙ yr⁻¹) are obtained. The SFH timescale also produces dependencies; for τ_SFH > 100 Myr (a larger value is closer to a constant SFH), we obtain M_* ∼ (1−10) × 10⁹ M_⊙ and SFR ∼ 10²−10³ M_⊙ yr⁻¹, and for τ_SFH < 100 Myr, M_* and SFR show larger variations. The metallicity is not constrained at all.

For HD2, the best-fit redshift is z = 12.3 (χ² = 4.2) and the 1σ range (Δχ² < 1) is 12.0 < z < 12.7. We find two types of high-redshift solutions. One is a very young starburst: an age less than 10 Myr, M_* ∼ 7 × 10⁹ M_⊙, SFR ∼ 10³−10⁴ M_⊙ yr⁻¹, and A _V ∼ 0.8. τ_SFH and metallicity are not constrained. The other case is a massive and relatively mature galaxy: an age greater than 100 Myr, M_* ∼ 1 × 10¹¹ M_⊙, SFR < 10² M_⊙ yr⁻¹, A _V < 0.5, and τ_SFH < 60 Myr. The metallicity is not constrained.

Although they are statistically less likely given the larger χ² values, the possible Balmer break solutions are as follows. For HD1, we obtain z ∼ 3.9, age of 0.3–1 Gyr, M_* ∼ (6−10) × 10⁹ M_⊙, SFR < 0.1 M_⊙ yr⁻¹, A _V < 0.5, τ_SFH < 0.1 Gyr, and Z > 0.004. For HD2, we obtain z ∼ 3.5, age of 0.4–0.7 Gyr, M_* ∼ 1 × 10¹⁰ M_⊙, SFR ∼ 0 M_⊙ yr⁻¹, A _V < 0.1, τ_SFH < 0.03 Gyr, and Z > 0.02. These stellar masses of ∼10¹⁰ M_⊙ are ∼10 times smaller than known passive galaxies at z ∼ 4 (Glazebrook et al. 2017; Tanaka et al. 2019; Valentino et al. 2020). Therefore, even these cases will be interesting to examine further spectroscopically in the future.

5.1. Selection Completeness

To derive the rest-frame UV luminosity function of the z ∼ 12–16 galaxies, we estimate the selection completeness by conducting Monte Carlo simulations. We first make mock SEDs of galaxies at 9.0 < z < 19.0 (steps of 0.1) with UV spectral slopes of−3.0 < β_UV < −1.0 (steps of 0.1). The intergalactic medium attenuation is taken into account by using a prescription of Inoue et al. (2014a), resulting in almost zero flux densities at a wavelength bluer than the Lyα break. We then calculate fluxes in each band by integrating mock SEDs through our 10 filters (grizyJHK _s (K)[3.6][4.5]) and scale to have apparent magnitudes of 23.0−25.0 mag in the K _s (K) band, whose central wavelength corresponds to ∼1500 Å at z ∼ 13. We then perturb the calculated fluxes by adding photometric scatters based on a Gaussian distribution with a standard deviation equal to the flux uncertainties in each band. We generate 1000 mock galaxies at each redshift with UV spectral slopes following a Gaussian distribution with a mean of β_UV = −2.0 and a scatter of Δβ_UV = 0.2 (Rogers et al. 2014; Bowler et al. 2020). Finally, we select z ∼ 12–16 galaxy candidates with the same color selection criteria and calculate the selection completeness as a function of the K _s (K)-band magnitude and redshift, C(m, z), averaged over the UV spectral slope. Figure 6 shows the calculated selection completeness in the COSMOS and SXDS fields. Our selection criteria can select sources at 12 ≲ z ≲ 16. The mean redshift from the simulation is z = 14.3 and 14.6 in the COSMOS and SXDS fields, respectively, but in this paper we adopt z = 12.8 (z ∼ 13), which is the average of the nominal redshifts for HD1 and HD2, as the mean redshift of our H-dropout sample. The selection completeness is ∼70% even for very bright (23.0 mag) galaxies because our color criteria are very strict in order to remove foreground interlopers (see Figure 1) and miss some intrinsically red (β_UV ≳ −1.8) z ≳ 12 galaxies.

Zoom In Zoom Out Reset image size

Based on the results of these selection completeness simulations, we estimate the survey volume per unit area as a function of the apparent magnitude (Steidel et al. 1999),

$$\begin{eqnarray}&&{V}_{\mathrm{eff}}(m)=\int C(m,z)\displaystyle \frac{{dV}(z)}{{dz}}{dz},\end{eqnarray} \tag{ 5 }$$

where C(m, z) is the selection completeness, i.e., the probability that a galaxy with an apparent magnitude m at redshift z is detected and satisfies the selection criteria, and dV(z)/dz is the differential comoving volume as a function of redshift.

5.2. Contamination

The space number density of the z ∼ 13 galaxies that are corrected for incompleteness and contamination is calculated with the following equation:

$$\begin{eqnarray}&&\psi (m)=\left[1-{f}_{\mathrm{cont}}\right]\displaystyle \frac{{n}_{\mathrm{raw}}(m)}{{V}_{\mathrm{eff}}(m)},\end{eqnarray} \tag{ 6 }$$

where n_raw(m) is the surface number density of selected galaxies in an apparent magnitude bin of m, and f_cont is a contamination fraction. We estimate the contamination fraction of foreground sources by conducting Monte Carlo simulations. As discussed in Section 4, the most likely contaminants are z ∼ 4 passive galaxies whose Balmer breaks mimic the Lyman break at z ∼ 13. Stellar contaminations are not expected to be dominant, given observed colors of stellar sources (Figure 1).

To investigate contamination from various z ∼ 4 passive galaxies, we prepare three types of mock SEDs at 3 ≤ z ≤ 5 based on (1) a classic spectrum of elliptical galaxies in Coleman et al. (1980), (2) model spectra with color distributions similar to real passive galaxies, and (3) the z ∼ 4 solutions from the SED fittings. In case 1, we use a spectrum of old elliptical galaxies in Coleman et al. (1980) as an input SED. In case 2, we first generate model spectra of galaxies by using PANHIT assuming a delayed-τ SFH, τ = 0.01, 0.03, 0.1, 0.3, and 1 Gyr; stellar age of 0.01−1.3 Gyr; metallicity of Z = 0.0001, 0.0004, 0.004, 0.008, 0.02, and 0.05; and dust attenuation of E(B−V) = 0−1 (steps of 0.05). Then, we calculate rest-frame near-UV−r and r−J colors of the models and compare these colors with those of observed passive galaxies in Davidzon et al. (2017). By selecting galaxies whose colors are consistent with the observed passive galaxies, we construct a set of passive galaxy SEDs that have realistic color distributions. In case 3, we use a spectrum of the z ∼ 4 passive galaxy solution in the SED fitting in Section 4. Since in this case we assume that all of the passive galaxies have the same SED as the z ∼ 4 solution, this case provides the most conservative estimate for the contamination fraction (i.e., the highest contamination fraction).

We then make mock SEDs redshifted to 3 ≤ z ≤ 5 (steps of 0.1) from the three types of the SEDs, calculate fluxes in each band, scale to have stellar masses of 10⁹ ≤ M_*/M_⊙ ≤ 10¹¹ (steps of 0.1 dex), and perturb the calculated fluxes by adding photometric scatters in the same manner as Section 5.1. We generate ∼1000 mock galaxies at each redshift and stellar mass bin and calculate the fraction of passive galaxies that satisfy our selection criteria in each bin. Finally, by integrating the product of the stellar mass function of passive galaxies in Davidzon et al. (2017) and the fraction of passive galaxies satisfying our selection criteria over the redshift and stellar mass, we calculate the number of passive galaxies at z ∼ 4 that are expected to be in our z ∼ 13 galaxy sample.

The expected numbers of passive galaxies in our sample are N_cont = 0.00, 0.12, and 1.36 in cases 1, 2, and 3, respectively. We estimate the contamination fraction f_cont by dividing the expected number of passive galaxies by the number of our z ∼ 13 candidates. The estimated contamination fractions are small in cases 1 and 2 (f_cont ∼ 0% and 6%, respectively) and f_cont ∼ 70% in case 3, where we assume that all of the passive galaxies have the same SED as the z ∼ 4 solution in the SED fitting as the extremely conservative case. Although the realistic simulation with the observed color distributions (i.e., case 2) indicates the very low contamination fraction, we adopt this very conservative estimate from case 3 for the UV luminosity function calculation. Note that even if we assume this conservative estimate as the prior, the z > 12 solutions for HD1 and HD2 in the SED fitting are still more likely than the z ∼ 4 solutions that give larger χ² values, as long as the true number density of z ∼ 13 galaxies is ≳10⁻⁸ Mpc⁻³ (comparable to our estimate in Section 5.3). On the other hand, if the true number density is ∼10⁻¹¹ Mpc⁻³ at z ∼ 13 (comparable to model predictions in Section 5.4), the z ∼ 4 solutions are more likely due to the higher number density of z ∼ 4 passive galaxies compared to that of z ∼ 13 galaxies.

5.3. UV Luminosity Function

We convert the number density of z ∼ 13 galaxies as a function of apparent magnitude, ψ(m), into the UV luminosity functions, Φ[M_UV(m)], which are the number densities of galaxies as a function of rest-frame UV absolute magnitude. We calculate the absolute UV magnitudes of galaxies from their apparent magnitudes in the K _s (K) band, whose central wavelength corresponds to ∼1500 Å at z ∼ 13, assuming a flat rest-frame UV continuum, i.e., constant f_ν , suggested by the SEDs of our galaxies. The 1σ uncertainty is calculated by taking into account the Poisson confidence limit (Gehrels 1986) on the expected number of galaxies at z ∼ 13 in our sample (N = 2 × (1−f_cont) ∼ 1).

Figure 7 shows the calculated UV luminosity function at z ∼ 13. The number density of our z ∼ 13 galaxies is $({3.7}_{-3.0}^{+8.4})\times {10}^{-8}\ {\mathrm{Mpc}}^{-3}\ {\mathrm{mag}}^{-1}$ at M_UV = −23.5 mag. This number density is comparable to that of bright galaxies at z ∼ 10 in Bowler et al. (2020), which is supported by the little evolution of the abundance of bright galaxies found by previous studies at z = 4–10 (Bowler et al. 2020; Harikane et al. 2022). Indeed, as shown in Figure 8, the number density of bright (M_UV < −23 mag) galaxies does not show significant redshift evolution from z ∼ 4 to z ∼ 13. In Figure 7, we also plot the number density of z ∼ 10 galaxies estimated from GN-z11, $({1.0}_{-0.8}^{+2.2})\times {10}^{-6}\ {\mathrm{Mpc}}^{-3}$. ¹⁷ These results and the spectroscopic confirmation of GN-z11 by Oesch et al. (2016) and Jiang et al. (2021) indicate that the bright end of the luminosity function at high redshift cannot be explained by the Schechter function with the exponential cutoff and is more consistent with the double-power-law function. Indeed, the number density of our z ∼ 13 galaxies is consistent with the double-power-law function with ${M}_{\mathrm{UV}}^{* }=-17.6$ mag, ϕ^* = 1.0 × 10⁻⁴ Mpc⁻³, α = −1.8, and β = −2.6 (Figure 7), which are derived by extrapolating the redshift evolution of the parameters in Bowler et al. (2020) to z = 13. ¹⁸ Although spectroscopic confirmation is needed, these results indicate that upcoming surveys will detect a number of galaxies at z > 10, which will be discussed in Section 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.4. Comparison with Models

Both theoretical and empirical models predict the UV luminosity function of galaxies at z > 10 (e.g., Dayal et al. 2014, 2019; Mason et al. 2015; Tacchella et al. 2018; Behroozi et al. 2019, 2020, Yung et al. 2019, 2020; see also Hutter et al. 2021). We compare the number densities at z ∼ 10 and 13 with predictions from these models in Figure 9. At z ∼ 10, the predictions roughly agree with the observed number densities for relatively faint galaxies (M_UV ≳ −21 mag), but the models underestimate the number densities of bright galaxies (M_UV ≲ −22 mag), albeit with large uncertainties in the observations. Similarly at z ∼ 13, the models cannot reproduce the observed number density of our z ∼ 13 galaxy candidates. These discrepancies indicate that the current models do not account for the rapid mass growth within the short physical time since the big bang.

Zoom In Zoom Out Reset image size

There are several possible physical processes to reconcile these discrepancies between the models and the observations at z ∼ 10–13. As discussed in Harikane et al. (2021b), less efficient mass quenching and/or lower dust obscuration than assumed in the models can explain the existence of these UV-bright galaxies. Active galactic nucleus (AGN) activity may also boost the UV luminosity in these galaxies. Previous studies indicate that the AGN fraction starts to increase at M_UV ≃ −22 mag (Ono et al. 2018; Stevans et al. 2018; Adams et al. 2020; Harikane et al. 2022; see also Piana et al. 2022). If we assume that the UV luminosities of HD1 and HD2 are solely powered by black holes, the inferred black hole masses are ∼10⁸ M_⊙, assuming accretion at the Eddington rate (Pacucci et al. 2022), in accordance with expectations for high-redshift quasars (see, e.g., Haiman & Menou 2000 and Willott et al. 2010). In addition, note that a ∼10⁸ M_⊙ black hole at z ∼ 12 could be the progenitor of z ∼ 7 quasars, as the growth time to reach a mass of 10⁹–10¹⁰ M_⊙, typical of z > 6 quasars detected thus far, is shorter than the cosmic time between z = 12–13 and z = 7, for an Eddington-limited accretion. It is also possible that the observed bright sources at z ∼ 10–13 are galaxies in a short-time starburst phase that is not captured in the models whose outputs are averaged over a time interval (see also Dayal et al. 2013). Finally, a top-heavier IMF would explain the discrepancies by producing more UV photons at the same stellar mass. It is possible that these bright galaxies (especially HD1 and HD2) are merging systems that are not resolved in the ground-based images. However, even if they are major mergers, the UV luminosity will decrease only by a factor of a few, which would not explain the discrepancy at z ∼ 13 (see also discussions in Harikane et al. 2022; Shibuya et al. 2022). In any case, if these bright z ∼ 10–13 galaxies are spectroscopically confirmed, the discrepancies will motivate the exploration of new physical processes that are responsible for driving the formation of these bright galaxies in the early universe.

5.5. Cosmic SFR Density

We calculate the cosmic SFR density at z ∼ 13 by integrating the rest-frame UV luminosity function. We use a double-power-law luminosity function at z = 13 with ${M}_{\mathrm{UV}}^{* }=-17.6$ mag, ϕ^* = 1.0 × 10⁻⁴ Mpc⁻³ α = −1.8, and β = −2.6, which is consistent with our number density measurement (see Section 5.3). We obtain the UV luminosity density by integrating the luminosity function down to −17 mag as in previous studies (e.g., Bouwens et al. 2015, 2020; Finkelstein et al. 2015; Oesch et al. 2018). We then convert the UV luminosity density to the SFR density by using the calibration used in Madau & Dickinson (2014) with the Salpeter (1955) IMF:

$$\begin{eqnarray}&&{\mathrm{SFR}}_{\mathrm{UV}}\ ({M}_{\odot }\,{\mathrm{yr}}^{-1})=1.15\times {10}^{-28}{L}_{\mathrm{UV}}\ (\mathrm{erg}\ {{\rm{s}}}^{-1}\ {\mathrm{Hz}}^{-1}).\end{eqnarray} \tag{ 7 }$$

This SFR density estimation is true only if the rest-frame UV galaxy identification is complete with respect to all the galaxy populations at z ∼ 13 (but see Fudamoto et al. 2020). The uncertainty of the SFR density is scaled from that of the number density measurement.

The estimated SFR density is ${\rho }_{\mathrm{SFR}}=({8.0}_{-6.6}^{+18.4})\times {10}^{-5}\ {M}_{\odot }{\mathrm{yr}}^{-1}\ {\mathrm{Mpc}}^{-3}$ at z = 12.8. We compare the SFR density with previous results in Figure 10. The estimated SFR density at z ∼ 13 is consistent with the fitting function in Harikane et al. (2021b), which is calibrated with recent observations at z > 4 showing a more rapid decline (∝10^−0.5(1+z)) than the extrapolation of the fitting function in Madau & Dickinson (2014) at z > 10. Furthermore, if we divide our sample into one galaxy at z ∼ 12 and another at z ∼ 13, given the low completeness at z = 12.3 (see Figure 6), the estimated SFR density would show a decrease from z ∼ 12 to 13, consistent with Harikane et al. (2021b). The estimate is also comparable to the edge of the range of the SFR density at z ≳ 13 expected from passive galaxies at z ∼ 6 in Mawatari et al. (2020a). Note that our estimated SFR density is dominated by relatively faint (M_UV ≳ −20 mag) galaxies (see Roberts-Borsani et al. 2021a for discussions at z ∼ 8–10). In this calculation, we need to assume the shape of the UV luminosity function, because there is no constraint on the number density of faint galaxies at z ∼ 13. Future JWST observations will measure the number density of faint z ∼ 13 galaxies and constrain the shape of the luminosity function combined with this work for bright galaxies, allowing a more robust measurement of the SFR density.

Zoom In Zoom Out Reset image size

There are several upcoming space-based missions that can search for z > 10 galaxies, such as JWST, Nancy Grace Roman Space Telescope (hereafter Roman), and Galaxy Reionization EXplorer and PLanetary Universe Spectrometer (GREX-PLUS). By taking advantage of the sensitivity of infrared observations in space, these missions are expected to detect galaxies at z > 10. In this section, we will discuss future prospects of these space missions based on our result of the z ∼ 12–16 galaxy search.

We here consider three missions, JWST, Roman, and GREX-PLUS, whose survey parameters are summarized in Table 4. Although Euclid has a remarkable capability to conduct wide-field surveys, it has a limited wavelength coverage to the H band and can only select sources up to z ∼ 10 (J-dropout). We detail the three missions and survey plans below.

Table 4. Expected Number of z ≳ 13 Galaxies Identified in Upcoming Surveys

Telescope	Survey	m5σ	Asurvey	N(z ∼ 13)	N(z ∼ 15)	N(z ∼ 17)
(1)	(2)	(3)	(4)	(5)	(6)	(7)
JWST	JADES/Deep	30.6/30.2	0.013	11–0.5	1−0	0−0
	JADES/Medium	29.7/29.3	0.053	15–0.5	1−0	0−0
	CEERS	29.0/29.2	0.027	3−0	0.5–0	0−0
	COSMOS-Web	−/28.1 a	0.6	19−0.5	2−0	0−0
	PRIMER	≤27.8/≤27.7 b	0.111	5−0	0.5−0	0−0
	NGDEEP	30.6/30.7	0.0027	2−0	1−0	0−0
	PANORAMIC	28.3/28.3	0.4	17−0.5	2−0	0−0
Roman	UltraDeep	29.6/28.6	1	250–7 c	5−0	⋯
	Deep	27.5/27.2	20	270–3 c	9−0	⋯
	HLS	26.7/-	2000	8441–33 c	⋯	⋯
GREX-PLUS	UltraDeep	27.7	1	18−0	1−0	0−0
	Deep	27.0	40	262−2	17−0	1−0
	Medium	26.0	200	300−0	17−0	1−0
	Wide	24.5	2000	322−0	14−0	0.5–0

Notes. Col. (1): telescope. Col. (2): planned survey. Col. (3): 5σ depth in the AB magnitude in the rest-frame UV band. In JWST, we quote depths of the F200W and F227W bands for the z ∼ 13 and z ∼ 15–17 galaxy selections, respectively. In Roman, depths of the F184 and K213 bands are quoted for the z ∼ 13 and z ∼ 15 galaxy selections, respectively. In GREX-PLUS, we quote depths in the F232 band for the z ∼ 13–17 galaxy selection. These depths are for point sources for simplicity, and the actual depths for high-redshift galaxies would be slightly shallower, except in PRIMER and NGDEEP, where resolved sources with typical sizes of high-redshift galaxies are assumed. Col. (4): survey area in deg². Cols. (5)–(7): expected number of galaxies at z ∼ 13, 15, and 17 identified in the survey assuming Δz = 1. Two values indicate the numbers in case A and case B (see text for details). Note that it may be difficult to select z ∼ 17 sources with Roman owing to the lack of observational bands redder than the K213 band. Because the current plan of the Roman HLS does not include the K213 imaging, we do not calculate the expected number of z ∼ 15 galaxies.

^a The COSMOS-Web survey will take NIRCam F115W, F150W, F277W, and F444W images. We use the depth of the F277W band for the z ∼ 13–17 galaxy selections. ^b PRIMER is a "wedding cake" survey that is composed of several surveys with different depths and areas (27.8/27.7 mag for $400\ {\mathrm{arcmin}}^{2}$, 28.3/28.3 mag for $300\ {\mathrm{arcmin}}^{2}$, and 28.8/28.8 mag for $33\ {\mathrm{arcmin}}^{2}$). ^c Although the expected number in the Roman HLS is larger than those in the other two Roman surveys, these three surveys will identify galaxies in different luminosities and are complementary to each other. Please see Figure 11 for luminosity ranges that each survey covers.

Download table as:  ASCIITypeset image

JWST is NASA's infrared space telescope that was launched on 2021 December 25. Thanks to its large 6.5 m diameter mirror, the sensitivity of JWST in the infrared is much higher than previous and current observational facilities. NIRCam is a near-infrared camera at 0.6–5.0 μm whose FOV is roughly $2\times 2.2\times 2.2\ {\mathrm{arcmin}}^{2}$ (Rieke et al. 2005). In this paper we consider the following six surveys using NIRCam: the JWST Advanced Deep Extragalactic Survey (JADES; Guaranteed Time Observation (GTO) program; Eisenstein et al. 2017a, 2017b), the Cosmic Evolution Early Release Science (CEERS) survey (ERS-1345; Finkelstein et al. 2017), the COSMOS-Web survey (GO-1727; Kartaltepe et al. 2021), Public Release IMaging for Extragalactic Research (PRIMER; GO-1837; Dunlop et al. 2021), the Next Generation Deep Extragalactic Exploratory Public (NGDEEP) survey (GO-2079; Finkelstein et al. 2021b), and the Parallel wide-Area Nircam Observations to Reveal And Measure the Invisible Cosmos (PANORAMIC) survey (GO-2514; Williams et al. 2021). These surveys will take deep NIRCam imaging data at ∼1–5 μm in ∼10–2000 arcmin² survey fields. Other high-redshift galaxy surveys are also planned in Cycle 1. For example, Through the Looking GLASS (ERS-1324; Treu et al. 2017) and the Ultra-deep NIRCam and NIRSpec Observations Before the Epoch of Reionization (UNCOVER; GO-2561; Labbe et al. 2021) will image the gravitational lensing cluster A2744 with NIRCam. The Webb Medium-Deep Field survey (GTO-1176; Windhorst et al. 2017) will use 110 hr to observe 13 medium-deep (28.4–29.4 mag) fields, including the JWST North Ecliptic Pole Time-domain Field (Jansen & Windhorst 2018).

Roman is NASA's optical to near-infrared space telescope whose launch is targeted around the mid-2020s. Although the size of the mirror is comparable to that of HST, Roman is expected to conduct wide-field surveys in the near-infrared by taking advantage of the wide FOV of its camera (0.28 deg²). The High Latitude Survey (HLS) will take images in the Y106, J129, H158, and F184 bands over the ∼2000 deg² sky, reaching ∼26.7 mag in F184. Two additional survey concepts are potentially possible: the Roman ultradeep field survey (hereafter UltraDeep; Koekemoer et al. 2019), which will take very deep images (reaching ∼30 mag in the bluer filters and 29.6 and 28.6 mag in F184 and K213, respectively) in a small area (≲1 deg²), and the Roman cosmic dawn survey (hereafter Deep; see also Rhoads et al. 2018), which will conduct a relatively wide and deep survey (∼20 deg², 27.5 and 27.2 mag in F184 and K213). ¹⁹ Both of these two possible surveys would cover a wavelength range up to 2.3 μm (the K213 filter), allowing us to select z ∼ 15 galaxies with the F184-dropout selection.

GREX-PLUS is a new 1.2 m class, cryogenic, wide-field infrared space telescope mission concept proposed to ISAS/JAXA for its launch around the mid-2030s. GREX-PLUS is planned to have a wide-field camera that will efficiently take wide and deep images at 2–10 μm, allowing us to select galaxies at z ∼ 10–17. Four surveys with different depths and areas (UltraDeep, Deep, Medium, and Wide) are now being planned.

For the surveys by these three missions, we calculate the expected number of detected galaxies at z ∼ 13, 15, and 17. For simplicity, we assume the survey volume of Δz = 1 and 100% completeness. We calculate the number of galaxies based on two cases of the rest-frame UV luminosity functions (case A and case B). One is the double-power-law function with a redshift evolution suggested in Bowler et al. (2020; case A), which is consistent with the number density of our z ∼ 13 galaxies. The other is the Schechter function with a density evolution, whose parameters are ${M}_{\mathrm{UV}}^{* }=-20.5$ mag, ϕ^* = 44.7 × 10^−0.6(1+z) Mpc⁻³, and α = −2.3 (case B). We assume a relatively rapid decrease of the ϕ^* parameter compared to the cosmic SFR density evolution at z > 10 (∝10^−0.5(1+z); Harikane et al. 2022), as a conservative estimate. These parameters are comparable to measurements for the z ∼ 10 UV luminosity function in Oesch et al. (2018). These two cases mostly cover the range of model predictions at z ∼ 13 as we show later in Figure 11. Given that these two cases are somewhat consistent with the model predictions and observations at z ∼ 10 and 13, we extrapolate these calculations to z ∼ 15 and 17 for reference. We integrate the luminosity functions down to the depths of detection bands presented in Table 4, and we calculate the number of detected galaxies at each redshift for each survey. Note that the 6σ depth is used for the z ∼ 15 galaxies identified with Roman to reduce the false-positive rate in the one-band (K213) detection, while the 5σ depth is adopted for the others because multiple bands can be used.

Zoom In Zoom Out Reset image size

Figures 11, 12, and 13 show the expected UV luminosity functions in case A and case B and the depth and volume of each survey at z ∼ 13, 15, and 17, respectively, with predictions from models at z ∼ 13 and 15. Table 4 summarizes the expected number of detected galaxies in each survey. JWST will identify galaxies up to z ∼ 15 in case A. In addition, JWST can conduct deep photometric and spectroscopic observations for relatively bright z > 10 galaxies identified in surveys with JWST and other telescopes, which will allow us to investigate physical properties (e.g., systemic redshift, stellar age, metallicity) in detail (e.g., Roberts-Borsani et al. 2021b). Roman will detect galaxies up to z ∼ 15 in case A and identify galaxies at z ∼ 13 even in case B thanks to the wide survey areas. GREX-PLUS may be able to push the redshift frontier to z ∼ 17 in case A. The wide-area surveys with Roman and GREX-PLUS can identify luminous z > 10 galaxies with ≲27 mag. These galaxies are bright enough to be followed up by spectroscopically with ALMA and JWST within a reasonable amount of observing time, to investigate the physical conditions of these galaxies in the early universe.

Zoom In Zoom Out Reset image size

In this paper we have presented our search for H-dropout LBGs at z ∼ 12–16. We have used the multiwavelength deep imaging data in the COSMOS and SXDS fields, including Subaru/HSC grizy, VISTA JHK _s , UKIRT JHK, and Spitzer/IRAC [3.6][4.5] images. Our major findings are summarized below:

• 1.  
  After the careful screening of foreground interlopers, we have identified two z ∼ 12–13 galaxy candidates, HD1 and HD2 (Figure 2). SEDs of these candidates show a sharp discontinuity between H and K _s (K) bands, nondetections in the grizyJ bands, and a flat continuum up to the [4.5] band, all of which are consistent with a z ∼ 12–13 galaxy. Photometric redshift analyses based on these SEDs indicate that the most likely redshifts are z > 12 for both HD1 and HD2.
• 2.  
  We calculate the number density of our galaxy candidates whose mean redshift is z ∼ 13 (Figure 7). The calculated number density at z ∼ 13 is comparable to that of bright galaxies at z ∼ 10 and consistent with the double-power-law luminosity function extrapolated to z = 13 assuming the redshift evolution in Bowler et al. (2020). These results support little evolution of the abundance of bright galaxies to z ∼ 13 as suggested by previous studies at z ∼ 4–10. Comparisons with theoretical and empirical models show that these models underestimate the number densities of bright galaxies at z ∼ 10–13, although the uncertainties in observations are large (Figure 9). The inferred cosmic SFR density is consistent with the rapid decrease at z > 10 with ∝ 10^−0.5(1+z) suggested by Harikane et al. (2021b) (Figure 10).
• 3.  
  We conducted ALMA follow-up observations targeting HD1. The obtained spectrum shows a ∼4σ tentative line-like feature at 237.8 GHz that is cospatial with the rest-frame UV emission, consistent with the O iii 88 μm emission line at z = 13.27 (Figure 5). Further spectroscopic efforts are needed to confirm the redshifts of HD1 and HD2.

Our results support the possibility that a number of bright galaxies exist at z > 10. If the UV luminosity function follows the double-power-law function consistent with the number density of the bright galaxies at z ∼ 10–13, upcoming space missions such as JWST, Roman, and GREX-PLUS will detect more than ∼10,000 galaxies at z ∼ 13–15 (Figures 11 and 12), and perhaps one to several at z ∼ 17 (Figure 13), allowing us to observe the first galaxy formation.

Zoom In Zoom Out Reset image size

We thank the anonymous referee for a careful reading and valuable comments that improved the clarity of the paper. We thank James Rhoads, Sangeeta Malhotra, Masami Ouchi, and the other members of the Roman cosmic dawn Science Investigation Team (SIT) for helpful discussions on the detectability of z > 10 galaxies with Roman. We are grateful to Caitlin Casey, James Dunlop, Steven Finkelstein, Christina Williams, and Rogier Windhorst for providing the expected depths in their JWST surveys, namely, COSMOS-Web, PRIMER, NGDEEP, PANORAMIC, and the Webb Medium-Deep Field survey, respectively. We thank Takashiro Morishita for bringing an error in Figure 6 in the earlier manuscript to our attention. This work was partly supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo, JSPS KAKENHI grant Nos. 17H06130, 19J01222, 20K22358, and 21K13953; the NAOJ ALMA Scientific Research Grant Codes 2018-09B and 2020-16B; and the black hole Initiative at Harvard University, which is funded by grants from the John Templeton Foundation and the Gordon and Betty Moore Foundation. T.H. was supported by Leading Initiative for Excellent Young Researchers, MEXT, Japan (HJH02007). P.D. and A.H. acknowledge support from the European Research Council's starting grant ERC StG-717001 ("DELPHI"). P.D. also acknowledges support from the NWO grant 016.VIDI.189.162 ("ODIN") and the European Commission's and University of Groningen's CO-FUND Rosalind Franklin program. A.Y. is supported by an appointment to the NASA Postdoctoral Program (NPP) at NASA Goddard Space Flight Center, administered by Oak Ridge Associated Universities under contract with NASA. F.P. acknowledges support from a Clay Fellowship administered by the Smithsonian Astrophysical Observatory.

This paper makes use of the following ALMA data: ADS/JAO. ALMA #2019.A.00015.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST, and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium.

The HSC collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office; the Ministry of Education, Culture, Sports, Science and Technology (MEXT); the Japan Society for the Promotion of Science (JSPS); Japan Science and Technology Agency (JST); the Toray Science Foundation; NAOJ; Kavli IPMU; KEK; ASIAA; and Princeton University.

This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org

This paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of the Center for Computational Astrophysics (CfCA), NAOJ.

The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, Johns Hopkins University, Durham University, the University of Edinburgh, the Queens University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Software: BEAGLE (Chevallard & Charlot 2016), PANHIT (Mawatari et al. 2020b), SExtractor (Bertin & Arnouts 1996), T-PHOT (Merlin et al. 2016).