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

We present two bright galaxy candidates at z ∼ 12 – 13 identi ﬁ ed in our H -dropout Lyman break selection with 2.3 deg 2 near-infrared deep imaging data. These galaxy candidates, selected after careful screening of foreground interlopers, have spectral energy distributions showing a sharp discontinuity around 1.7 μ m, a ﬂ at continuum at 2 – 5 μ m, and nondetections at < 1.2 μ m in the available photometric data sets, all of which are consistent with a z > 12 galaxy. An ALMA program targeting one of the candidates shows a tentative 4 σ [ O III ] 88 μ m line at z = 13.27, in agreement with its photometric redshift estimate. The number density of the z ∼ 12 – 13 candidates is comparable to that of bright z ∼ 10 galaxies and is consistent with a recently proposed double-power-law luminosity function rather than the Schechter function, indicating little evolution in the abundance of bright galaxies from z ∼ 4 to 13. Comparisons with theoretical models show that the models cannot reproduce the bright end of rest-frame ultraviolet luminosity functions at z ∼ 10 – 13. Combined with recent studies reporting similarly bright galaxies at z ∼ 9 – 11 and mature stellar populations at z ∼ 6 – 9, our results indicate the existence of a number of star-forming galaxies at z > 10, which will be detected with upcoming space missions such as the James Webb Space Telescope, Nancy Grace Roman Space Telescope, and GREX-PLUS.


INTRODUCTION
hari@icrr.u-tokyo.ac.jpObserving the first galaxy formation is one of the main goals in the modern astronomy.One of the most straightforward approaches to achieve this goal is to observe forming galaxies directly in the early universe.arXiv:2112.09141v4[astro-ph.GA] 12 Feb 2022 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 serious problem to form blackholes as massive as ∼ 10 9 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 is not only the simplest frontier of the knowledge of human beings but also has a 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 2 , several studies using HST report very luminous Lyman break galaxies (LBGs) at z ∼ 9 − 10 more frequently than the expectation from a Schechtershape 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 a few square-degree near-infrared imaging surveys with Visible and Infrared Survey Telescope for Astronomy (VISTA) and UK Infrared Telescope (UKIRT) such as UltraVISTA (McCracken et al. 2012), the UKIRT In-fraRed 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 than a standard Schechter function (Stefanon et al. 2017(Stefanon et al. , 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. 2021b), implying little evolution of the number density of bright galaxies at z ∼ 4 − 10 (Bowler et al. 2020;Harikane et al. 2021b).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 age is ∼ 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. 2020b).
Motivated by these recent works, we search for Hband dropout (H-dropout) LBGs whose plausible redshifts are z ∼ 12 − 16. 1 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 datasets 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 James Webb Space Telescope (JWST).
This paper is organized as follows.We describe photometric datasets 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 summarize our findings in Section 7. Throughout this paper, we use the Planck cosmological parameter sets of the TT, TE, EE+lowP+lensing+ext result (Planck Collaboration et al. 2016): Ω m = 0.3089, Ω Λ = 0.6911, Ω b = 0.049, h = 0.6774, and σ 8 = 0.8159.All magnitudes are in the AB system (Oke & Gunn 1983).

Dataset
We use deep and wide photometric datasets available in the COSMOS (Scoville et al. 2007) and SXDS (Furu-

Selection of H-dropout Galaxies
We construct multi-band 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 (∼ 1. 7), 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. (2018Harikane et al. ( , 2019)).As high-resolution prior images in the T-PHOT run, we use HSC grizy stacked images whose PSF is ∼ 0. 7. 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 multi-band photometric catalogs constructed above.We adopt the following color selection criteria in the COSMOS and SXDS fields, respectively: COSMOS: SXDS: 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.
To remove foreground interlopers further, we conduct a photometric redshift analysis using BEAGLE (Chevallard & Charlot 2016).We adopt a constant star forma-0.01.0

Galactic stars
Figure 1.Two-color diagrams to select H-dropout galaxies.The left and right panels show two-color diagrams in the COSMOS and SXDS fields, respectively.The red lines indicate color criteria that we use to select H-dropout galaxies (Equations ( 1)-( 4)), and the red squares are the selected candidates, HD1 and HD2.The black solid lines are colors of star-forming galaxies at z ≥ 12 calculated with BEAGLE (Chevallard & Charlot 2016) with τV = 0.0 and 0.4 (corresponding to UV spectral slopes of βUV −2.4 and −1.9, respectively) as a function of redshift.The circles on the line show their redshifts with an interval of ∆z = 0.1.The blue circles are z = 0 − 7 sources spectroscopically identified (Laigle et al. 2016;Mehta et al. 2018).The dotted, dashed, and dot-dashed lines are, respectively, typical spectra of elliptical, Sbc, and irregular galaxies (Coleman et al. 1980) redshifted from z = 0 to z = 7.The black stars indicate Galactic dwarf stars taken from Patten et al. (2006) and Kirkpatrick et al. (2011).tion history with the Chabrier (2003) initial mass function (IMF), stellar ages of 10 6 , 10 7 , and 10 8 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 ∆χ 2 > 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.
HD1 is also found in the COSMOS2020 catalog (Weaver et al. 2021).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 re-measure magnitudes by using a larger aperture in the original IRAC images before the T-PHOT run, the magnitudes become brighter due to a neighboring source located ∼ 3. 5 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. 2021a).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.

ALMA FOLLOW-UP OBSERVATION
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 [Oiii]88µm emission lines in high redshift galaxies (e.g., Inoue et al. 2016, Laporte et al. 2017, 2021, Carniani et al. 2017, 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 [Oiii]88µm line using four tuning setups with the 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 typical beam size is ∼ 0. 4 − 0. 7. 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 −1 .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 1. 0-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 mo- ment 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 3 absorption.If this feature is the [Oiii]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 −1 in a full width at half maximum; FWHM) is in fact comparable to similarly bright LBGs at z ∼ 6 (Harikane et al. 2020).The emission feature is co-spatial with the rest-frame UV emission in the K s -band image (the bottom panel of Figure 5).The integrated line flux is 0.24 ± 0.06 Jy km s −1 or (1.9 ± 0.5) × 10 −18 erg s −1 cm −2 , and the line luminosity is The line luminosity is very small compared to the UV luminosity.Since the UV luminosity of HD1 is This ratio is the smallest among the galaxies observed in [Oiii]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 [Oiii]-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 followup 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. 2021a), which allows us to examine a wider redshift range than ALMA.
The dust continuum of HD1 remains non-detection, 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 −1 .Assuming that HD1 is not resolved in this observation, we obtain the 3σ upper limit on the dust continuum of < 24 µJy.

SED FITTING
To examine the photometric redshifts of HD1 and HD2 more carefully, we perform a comprehensive SED fitting analysis from optical to sub-millimeter (sub-mm) wavelength using PANHIT (Mawatari et al. 2020a).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 sub-mm range into account.PANHIT deals with the upper limits for non-detection 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 sub-mm 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 star formation history (SFH) covering a wide range of histories including a short time-scale 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 black-body 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.
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 χ 2 values than the z > 12 solutions.These are Balmer break galaxy so- lutions, and the dust temperature does not affect them because these solutions have very weak or no dust emission.Another type of possible solutions is dusty Hα emitters at z ∼ 2, although these solutions are not supported by the non-detections in the far-infrared and sub-mm 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 χ 2 than those of the Balmer break ones at z ∼ 4, while the lower, more normal dust temperature cases do not favor this type of the solution.Moreover, even the 80 K case is significantly less likely compared to the solutions at z > 12 (∆χ 2 > 4).For HD1, the best-fit redshift is z = 15.2 (χ 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 χ 2 = 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 9 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 9 M and SF R ∼ 10 2−3 M yr −1 , respectively.For an age of 10-100 Myr (> 100 Myr), M * ∼ (1 − 10) × 10 9 M and SF R ∼ 10 2 M yr −1 (M * ∼ (10 − 100) × 10 9 M and SF R < 10 2 M yr −1 ) 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 9 M and SF R ∼ 10 2−3 M yr −1 , and for τ SFH < 100 Myr, M * and SF R show larger variations.The metallicity is not constrained at all.

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 IGM attenuation is taken into account by using a prescription of Inoue et al. (2014a), resulting almost zero flux densities at the 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.0mag 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) that 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.
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), 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.

Contamination
The space number density of the z ∼ 13 galaxies that are corrected for incompleteness and contamination is calculated with the following equation: 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-τ star formation history, τ = 0.01, 0.03, 0.1, 0.3, and 1 Gyr, the 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 N U V − 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 9 ≤ M * /M ≤ 10 11 (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, 1.36 in cases of 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 of 1 and 2 (f cont ∼ 0% and 6%, respectively), and f cont ∼ 70% in case 3, where we assume 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 χ2 values, as long as the true number density of z ∼ 13 galaxies is 10 −8 Mpc −3 (comparable with our estimate in Section 5.3).On the other hand, if the true number density is ∼ 10 −11 Mpc −3 at z ∼ 13 (comparable with 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.

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 is 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 the calculated UV luminosity function at z ∼ 13.The number density of our z ∼ 13 galaxies is (3.7 +8.4  −3.0 ) × 10 −8 Mpc −3 mag −1 at M UV = −23.5 mag.This number density is comparable to that of bright Figure 7. Rest-frame UV luminosity functions at z ∼ 13 and z ∼ 10.The red circle shows the number density of our z ∼ 13 galaxy candidates.The black symbols and the gray shaded region are measurements at z ∼ 10 from the literature (diamond: McLeod et al. 2016, square: Oesch et al. 2018, pentagon: Morishita et al. 2018, circle: Bowler et al. 2020, shade: Finkelstein et al. 2021a).The green star is the number density of GN-z11 (see text).Note that the data point of Bowler et al. (2020) (GN-z11) is horizontally (vertically) offset by −0.2 mag (+0.03 dex) for clarity.The gray dashed line is the Schechter function fit (Bouwens et al. 2016), whereas the gray and red solid lines are the double power-law functions at z ∼ 10 and 13, respectively, whose parameters are determined by the extrapolation from lower redshifts in Bowler et al. (2020).
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. 2021b).Indeed as shown in Figure 8, the number density of bright (M UV < −23 mag) galaxies do 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 +2.2 −0.8 ) × 10 −6 Mpc −3 . 2 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 * UV = −17.6 mag, φ * = 1.0 × 10 −4 Mpc −3 , α = −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. 3 Although spectroscopic confirmation is needed, these results indicate that upcoming surveys will detect a number of galaxies at z > 10, which will be discussed later in Section 6.

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,   3 The extrapolated double power-law luminosity function from Bowler et al. (2020) predicts a ∼ 2 times higher number density at z ∼ 11 − 12 than those at z ∼ 10 and 13 in the magnitude regime of M UV −23.5 mag, but still consistent with our z ∼ 13 estimate within the errors.The detail of such a redshift evolution is beyond the scope of this paper.
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.
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.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. 2021b, see also Piana et al. 2021).If we assume that the UV luminosities of HD1 and HD2 are solely powered by black holes, the inferred black hole masses are ∼ 10 8 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 andWillott et al. 2010).In addition, note that a ∼ 10 8 M black hole at z ∼ 12 could be the progenitor of z ∼ 7 quasars, as the growth time to reach a mass of 10 9−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 source 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. 2021b andShibuya et al. 2021).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.Figure 9.Comparison with predictions from theoretical and empirical models at z ∼ 10 (left) and z ∼ 13 (right).The black symbols and the gray shaded region are measurements at z ∼ 10 from the literature (symbols are the same as in Figure 7) and at z ∼ 13 from this study.The blue lines show predictions from models (solid: Dayal et al. 2014, 2019, dotted: Yung et al. 2019, 2020, dot-dashed: Behroozi et al. 2019, 2020, dashed: Mason et al. 2015).
(7) This SFR density estimation is true only if the restframe UV galaxy identification is complete with respect to the all 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 ρ SFR = (8.0+18.4 −6.6 ) × 10 −5 M yr −1 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. (2020b).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, be-cause 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.

FUTURE PROSPECT
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.
JWST is NASA's infrared space telescope that was launched on 2021 December 25th.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 × 2.2 × 2.2 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,b), the Cosmic Evolution Early Release Science (CEERS) survey (ERS-1345;Finkelstein et al. 2017) (GTO-1176;Windhorst et al. 2017) will use 110 hours to observe 13 medium-deep (28.4 − 29.4 mag) fields, including the James Webb Space Telescope North Ecliptic Pole Timedomain Field (Jansen & Windhorst 2018).
Roman is NASA's optical to near-infrared space telescope whose launch is targeted around 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 2 ).The high latitude survey (HLS) will take images in the Y106, J129, H158, and F184bands over the ∼ 2000 deg 2 sky reaching ∼ 26.7 mag in F184.Two additional survey concepts are potentially possible; the Roman ultra-deep 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 2 ), 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 2 , 27.5 and 27.2 mag in F184 and K213). 4oth 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, widefield infrared space telescope mission concept proposed to ISAS/JAXA for its launch around 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 (Ultra-Deep, 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).The 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 * UV = −20.5 mag, φ * = 44.7 × 10 −0.6(1+z) Mpc −3 , and α = −2.3(case B). (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.(4) Survey area in deg 2 .( 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 due 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.* 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.
§ Although the expected number in the Roman HLS is larger than 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.
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. 2021b), 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 pre-sented in Table 4, and 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.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 m UV Figure 11.Predicted UV luminosity functions and depths and volumes of upcoming surveys.The red solid and dashed lines are the double power-law and Schechter functions (case A and case B in the text), respectively, whose parameters are determined by the extrapolation from lower redshifts (see text for details).The black, blue (left) and green (right) lines indicate expected depths and volumes of upcoming surveys with JWST, Roman, and GREX-PLUS, respectively.The gray lines show predictions from models (solid: Dayal et al. 2014, 2019, dotted: Yung et al. 2019, 2020, dot-dashed: Behroozi et al. 2019, 2020, dashed: Mason et al. 2015, same as the right panel of Figure 9).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 widearea 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.didates, HD1 and HD2 (Figure 2).SEDs of these candidates show a sharp discontinuity between H and K s (K)-bands, non-detections in the grizyJbands, 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 al. (2020).results support little evolution of the abundance of galaxies to z ∼ suggested 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 [Oiii]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 ∼ 10000 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.
(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 programme 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

Figure 2 .Figure 3 .
Figure2.Left: Optical-to-near-infrared SEDs of our H-dropout galaxy candidates, HD1 (top) and HD2 (bottom).The red symbols with error-bars are measured flux densities or the 2σ upper limits.The blue curve shows the best-fit model of an LBG at z > 12 in the SED fitting, and the gray curve shows a passive galaxy solution at z ∼ 4 (see Section 4).The upper panels show 10 × 10 images.The "[3.6] resid" and "[4.5] resid" are residual images after subtracting nearby objects with T-PHOT(Merlin et al. 2016).Middle: The same as the left panel but for the far-infrared-to-submm range.The curve shows the modified black-body function with the temperature of 50 K and the emissivity index of β dust = 2.0.Right: χ 2 value as a function of the redshift.The best-fit models are found at z > 12.

Figure 5 .
Figure 5. Top: ALMA spectrum showing the 4σ line-like feature at 237.8 GHz.This hints for the [Oiii]88µm line at z = 13.27.Bottom: Integrated intensity of the 4σ feature in HD1 overlaid on the VISTA Ks band image.This moment 0 map is made with the CASA task immoments, by integrating over 700 km s −1 covering most of the velocity range of the line emission (> 1.5×FWHM).The solid (dotted) lines show +1.5, +2.5, and +3.5σ (−1.5, −2.5, and −3.5σ) contours.The emission is co-spatial with the rest-frame UV emission in the Ks-band image.

Figure 6 .
Figure 6.Completeness estimated in our Monte Carlo simulations.The top and bottom panels are results for the COS-MOS and SXDS fields, respectively.The red, orange, yellow, green, blue, and purple lines show the completeness for Ks(K) = 23.0,23.4,23.8, 24.2, 24.6, and 25.0 mag sources, respectively.

Figure 8 .
Figure 8. Evolution of the rest-frame UV luminosity functions from z ∼ 4 to z ∼ 13.The red circle shows the number density of our z ∼ 13 galaxy candidates, and the grey, brown, purple, blue, green, yellow, and orange symbols show results at z ∼ 4, 5, 6, 7, 8, 9, and 10, respectively.The circles at z ∼ 4 − 7 are galaxy number densities from Harikane et al. (2021b), and those at z ∼ 8 − 10 are taken from Bowler et al. (2020).The squares show results taken from Bouwens et al. (2021) and Oesch et al. (2018) at z ∼ 4 − 9 and z ∼ 10, respectively.The diamond is a result in McLeod et al. (2016).The lines show double power-law functions in Harikane et al. (2021b) at z ∼ 4 − 7 and Bowler et al. (2020) at z ∼ 8 − 13.Note that the data point of Bowler et al. (2020) at z ∼ 10 is horizontally offset by −0.2 mag for clarity.
Figure 10.Evolution of the cosmic SFR density.The red circle is our result at z ∼ 13 estimated by integrating the double power-law luminosity function down to −17 mag.The black circles are observed cosmic SFR densities taken from Madau & Dickinson (2014), Finkelstein et al. (2015), McLeod et al. (2016), and Bouwens et al. (2020).The blue and gray dashed curves represent the fits in Harikane et al. (2021b, their Equation (60)) and Madau & Dickinson (2014, extrapolated at z > 8), respectively.All results are converted to use the Salpeter (1955) IMF (Equation (7)).
, 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 2 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 gravita-tional lensing cluster Abell 2744 with NIRCam.The Webb Medium-Deep Field survey

Figure 12 .
Figure 12.Same as Figure 11 but for z ∼ 15.The 6σ depth is used for Roman to reduce the false positive rate in the one-band (K213) detection.

Figure 13 .
Figure 13.Same as Figure 11 but for z ∼ 17.We only plot the survey depths and volumes of JWST and GREX-PLUS, which cover the wavelength redder than 2.3 µm (the Lyα break at z ∼ 17).

Table 1 .
5σ Limiting Magnitude of Imaging Data Used in This Study Note-R.A. and decl.are the central coordinates of the survey field.The 5σ limiting magnitudes are measured in 1. 5, 2. 0, and 3. 0-diameter apertures in grizy, JHKs(K), and [3.6][4.5]images, respectively, taken from Harikane et al. (2021b), release notes of UltraVISTA DR4 a (McCracken et al. 2012) and UKIDSS DR11 b (Lawrence et al. 2007), and Harikane et al. (2018).The values for the JHKs bands in the UD-COSMOS field represent limiting magnitudes in the ultra-deep and deep stripes.
(Lawrence et al. 2007)hara et al. , 2019) ) of the deep near-infrared data.Specifically, we use optical grizy images obtained in the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) Survey(Aihara et al. 2018(Aihara et al. , 2019) )public data release 2 (PDR2), near-infrared JHK s /K images of the UltraVISTA DR4(McCracken et al. 2012) and UKIDSS UDS DR11(Lawrence et al. 2007)in the COSMOS and SXDS fields, respectively, and Spitzer/IRAC [3.6] and [4.5] images obtained in the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH).Typical 5σ limiting magnitudes of these imaging data are presented in Table1.Since the COS-MOS field has the ultra-deep and deep stripes with different depths in the near-infrared images, we use the limiting magnitude of each stripe depending on the location of the source of interest.

Table 2 .
Photometry of Our H-Dropout Galaxy Candidates RedshiftFigure4.Full ALMA spectrum of HD1.This spectrum is extracted from a 1. 0-radius circular aperture centered on the coordinate of HD1.No obvious emission line is identified at > 5σ, but there is a 4σ line-like feature at 237.8 GHz (the red arrow), where no severe atmospheric O3 absorption exists (the gray shades).

Table 3 .
Physical Properties of Our H-Dropout Candidates

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