We conduct a comprehensive study on dropout galaxy candidates at z  916 using the first 90 arcminSUP2SUP James Webb Space Telescope JWST Near Infrared Camera images taken by the early release observations ERO and early release science programs. With the JWST simulation images we find that a number of foreground interlopers are selected with a weak photoz determination  SUP2SUP gt 4. We thus carefully apply a secure photoz selection criterion  SUP2SUP gt 9 and conventional color criteria with confirmations of the ERO Near Infrared Spectrograph spectroscopic redshifts and obtain a total of 23 dropout galaxies at z  916 including two candidates at zmathrmphot16.250.460.24 and 16.410.550.66 . We perform thorough comparisons of dropout galaxies found in our work with recent JWST studies and conclude that our galaxy sample is reliable enough for statistical analyses. We derive the UV luminosity functions at z  916 and confirm that our UV luminosity functions at z  9 and 12 agree with those determined by other Hubble Space Telescope and JWST studies. The cosmic star formation rate SFR density decreases from z  9 to 12 and perhaps to 16 but the densities at z  1216 are higher than the constant star formation efficiency model. Interestingly there are six bright galaxy candidates at z  1016 with M SUBUVSUB lt 19.5 mag and M SUBSUB  10SUP89SUP M SUBSUB. Because a majority 80 of these galaxies show no signatures of active galactic nuclei in their morphologies the high cosmic SFR densities and the existence of these UVluminous galaxies are explained by the lack of suppression of star formation by the UV background radiation at the prereionization epoch andor an efficient UV radiation production by a topheavy initial mass function with Population IIIlike star formation.
['Galaxy evolution', 'Galaxy formation', 'High-redshift galaxies', '594', '595', '734', 'Astrophysics - Astrophysics of Galaxies']
A Comprehensive Study of Galaxies at z   916 Found in the Early JWST Data Ultraviolet Luminosity Functions and Cosmic Star Formation History at the Prereionization Epoch
{'A Comprehensive Study on Galaxies at z ∼ 9 -16 Found in the Early JWST Data: UV Luminosity Functions and Cosmic Star-Formation History at the Pre-Reionization Epoch': 'Yuichi Harikane, 1 Masami Ouchi, 2, 1, 3 Masamune Oguri, 4, 5 Yoshiaki Ono, 1 Kimihiko Nakajima, 2 Yuki Isobe, 1, 6 Hiroya Umeda, 1, 6 Ken Mawatari, 2 and Yechi Zhang 1, 7  \n1 Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan 2 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan 3 Kavli Institute for the Physics and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba 277-8583, Japan 4 Center for Frontier Science, Chiba University, 1-33 Yayoi-cho, Inage-ku, Chiba 263-8522, Japan 5 Department of Physics, Graduate School of Science, Chiba University, 1-33 Yayoi-Cho, Inage-Ku, Chiba 263-8522, Japan 6 Department of Physics, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  \n- 7 Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-0033, Japan  \n(Accepted for Publication in ApJS)', 'ABSTRACT': 'We conduct a comprehensive study on dropout galaxy candidates at z ∼ 9 -16 using the first 90 arcmin 2 JWST/NIRCam images taken by the early release observations (ERO) and early release science (ERS) programs. With the JWST simulation images, we find that a number of foreground interlopers are selected with a weak photoz determination (∆ χ 2 > 4). We thus carefully apply a secure photoz selection criterion (∆ χ 2 > 9) and conventional color criteria with confirmations of the ERO NIRSpec spectroscopic redshifts, and obtain a total of 23 dropout galaxies at z ∼ 9 -16, including two candidates at z phot = 16 . 25 +0 . 24 -0 . 46 and 16 . 41 +0 . 66 -0 . 55 . We perform thorough comparisons of dropout galaxies found in our work with recent JWST studies, and conclude that our galaxy sample is reliable enough for statistical analyses. We derive the UV luminosity functions at z ∼ 9 -16, and confirm that our UV luminosity functions at z ∼ 9 and 12 agree with those determined by other HST and JWST studies. The cosmic star-formation rate density decreases from z ∼ 9 to 12, and perhaps to 16, but the densities at z ∼ 12 -16 are higher than the constant star formation efficiency model. Interestingly, there are six bright galaxy candidates at z ∼ 10 -16 with M UV < -19 . 5 mag and M ∗ ∼ 10 8 -9 M glyph[circledot] . Because a majority ( ∼ 80%) of these galaxies show no signatures of AGNs in their morphologies, the high cosmic star-formation rate densities and the existence of these UV-luminous galaxies are explained by no suppression of star-formation by the UV background radiation at the pre-reionization epoch and/or an efficient UV radiation production by a top-heavy IMF with Population III-like star formation.  \nKeywords: galaxies: formation - galaxies: evolution - galaxies: high-redshift', '1. INTRODUCTION': "One of the most important goals in astronomy today is to understand galaxy formation from their birth stage to current stage (Stark 2016; Dayal & Ferrara 2018; Ouchi et al. 2020; Robertson 2021). To accomplish the goal, observations for present galaxies to first galaxies are key to revealing the entire process of galaxy formation, while observations of early high redshift galaxies, especially  \nhari@icrr.u-tokyo.ac.jp  \nfirst galaxies, are missing (e.g., Zackrisson et al. 2011; Nakajima & Maiolino 2022).  \nOver the past 2-3 decades, large telescopes have driven observational studies of galaxy formation with millions of galaxies at a redshift up to z ∼ 10 since the start of deep imaging observations represented by the legendary Hubble Deep Field project with Hubble Space Telescope (HST; Williams et al. 1996). To date, deepfield imaging observations have reached detection limits of glyph[similarequal] 30 mag in the wavelength range of 0 . 4 -1 . 6 µ m with HST/ACS and WFC3 instruments in the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006, see  \nBouwens et al. 2021 and references therein) with the moderately deep ultra-violet (UV) extension, UVUDF (0 . 2 -0 . 4 µ m; Windhorst et al. 2011; Teplitz et al. 2013). Albeit with shallower detection limits of ∼ 26 -29 mag, HST GOODS, COSMOS, and CANDELS, and the associated parallel-field programs have covered a total area of square degrees in the blank fields (Giavalisco et al. 2004; Scoville et al. 2007; Grogin et al. 2011; Koekemoer et al. 2011). Complementary ground-based observations of Subaru Hyper Suprime-Cam survey have been completed optical imaging covering ∼ 1000 deg 2 with ∼ 26 mag depth (Aihara et al. 2022; see also Harikane et al. 2022b), while the ground-based nearinfrared (1 -2 µ m) and Spitzer Space Telescope imaging (3 -8 µ m) are limited to a total of few deg 2 with the similar or moderate depths of ∼ 25 -26 mag. Beyond deep imaging in blank fields, the HST programs, Hubble Frontier Fields (HFF) and Reionization Lensing Cluster Survey (RELICS), target 6 and 41 massive clusters, respectively, with depths of ∼ 26 -29 mag to study faint high redshift galaxies behind the clusters, exploiting gravitational lensing magnification (Lotz et al. 2017; Coe et al. 2019). These deep imaging data provide more than 4 million photometrically-selected dropout galaxies at z ∼ 4 -10 (Bouwens et al. 2021; Harikane et al. 2022b) and up to z ∼ 13 (Harikane et al. 2022a). Albeit with much small high redshift galaxy samples, spectroscopic observations confirm galaxies up to z = 9 . 1 with ALMA (Hashimoto et al. 2018) and z = 11 . 0 by HST/WFC3 grism and Keck/MOSFIRE spectroscopy (Oesch et al. 2016; Jiang et al. 2021). Star formation in even higher redshift ( z glyph[greaterorsimilar] 14) is discussed based on Balmer break galaxy candidates at z ∼ 6 (Mawatari et al. 2020b).  \nWith the galaxy samples photometrically selected in the rest-frame UV wavelengths, a number of studies have derived rest-frame UV luminosity functions reaching up to z ∼ 10 -13. The UV luminosity functions show the redshift evolution from z ∼ 3 to 10 with a decrease of the normalization φ ∗ and an increase of the faint-end slope α , and no evolution of characteristic luminosity L ∗ on the basis of the Schechter function parameterization (Bouwens et al. 2015, 2021; Finkelstein et al. 2015b; Ishigaki et al. 2018). At z ∼ 4 and above, there are claims that the bright-end of the UV luminosity function is explained with the double power-law function, but not with the Schechter function, due to the excessive number of bright galaxies (Bowler et al. 2014, 2020; Ono et al. 2018; Stevans et al. 2018; Harikane et al. 2022b). Such bright galaxy population includes the galaxies with spectroscopic redshifts at z ∼ 10 (Oesch et al. 2016; Jiang et al. 2021) and perhaps galaxy candidates at z ∼ 13 (Harikane et al. 2022a), while it is not clearly concluded  \nwith the sufficient statistical accuracy and the spectroscopic confirmations (e.g., Stefanon et al. 2019; RojasRuiz et al. 2020; Bowler et al. 2020).  \nOver the cosmological volumes, the redshift evolution of cosmic star formation rate (SFR) density is revealed with the UV luminosity function measurements, and shows a monotonic decrease from z ∼ 3 to z ∼ 10 with a small contribution of dusty starbursts at z glyph[greaterorsimilar] 6 (e.g., Bouwens et al. 2022a; Barrufet et al. 2022). The UV luminosity function measurements provide the physical picture of galaxy formation over the redshift range of z ∼ 0 -10, tying galaxies and dark-matter halos via abundance matching techniques (e.g. Behroozi et al. 2013, 2019; Moster et al. 2013, 2018; Finkelstein et al. 2015a). There is an increasing trend of stellar-to-halo mass ratio towards high redshift for a given halo mass (e.g., Behroozi et al. 2013, 2019; Harikane et al. 2016, 2018), which is consistent with the original idea of the galaxy-formation downsizing picture (Cowie et al. 1988). The galaxy-dark matter halo connection probed by the clustering analysis indicates that the star formation efficiency, defined by the ratio of the SFR to the dark matter accretion rate, SFR/ ˙ M h , is almost constant across the redshift of z ∼ 2 -7 given the dark matter halo mass (Harikane et al. 2018, 2022b), and the constant star formation efficiency model can reproduce the evolutional trend of the cosmic SFR density (e.g., Bouch'e et al. 2010, Mason et al. 2015, Harikane et al. 2018, 2022b, Tacchella et al. 2018).  \nThe UV luminosity function measurements, especially at the faint end, are clue to understanding galaxy formation (Yue et al. 2016) as well as cosmic reionization (Robertson 2021), where abundant faint star-forming galaxies are thought to be sources of cosmic reionization. The faint-end ( glyph[greaterorsimilar] -15 mag) UV luminosity function at z ∼ 6 -10 is probed with galaxies behind massive clusters such with the HFF data via gravitational lensing magnification (Atek et al. 2015, 2018; Ishigaki et al. 2015, 2018; Livermore et al. 2017; Laporte et al. 2016; Oesch et al. 2018), while the faint-end slopes and luminosity function turnovers are poorly constrained, due to the limited statistics and lensing magnification systematics (Bouwens et al. 2017, 2022a; Kawamata et al. 2018; Yang et al. 2022).  \nHere the James Webb Space Telescope (JWST) was launched in the end of 2021, and just started its operation in the early 2022. The first data sets of JWST were released on July 12, 2022, taken by the early release observations (ERO) whose targets include a massive cluster SMACS J0723.3-7327 (SMACS J0723, z = 0 . 39) and Stephan's Quintet. The ERO imaging data taken with NIRCam (Rieke et al. 2005) are deep enough to  \ndetect high redshift galaxies with the depths of ∼ 30 mag, and multi-band data covering glyph[greaterorsimilar] 2 µ m wavelengths allow us to detect galaxies at the previously unreachable redshift range up to z ∼ 20. Rest-frame optical emission at z glyph[greaterorsimilar] 10 is redshifted to the mid infrared bands and can be covered with the Mid-Infrared Instrument (MIRI; Bouchet et al. 2015). The ERO spectroscopic data of NIRSpec (Jakobsen et al. 2022) taken in the multi-object spectroscopy mode confirmed galaxies up to z = 8 . 5 with rest-frame optical lines in the 2 -5 µ m wavelengths. The slit-less spectroscopy of NIRISS (Doyon et al. 2012) supplement spectroscopic redshift determinations in the wavelength range of ∼ 1 -2 µ m. All of these data sets are revolutionizing the galaxy formation studies. The JWST observatory subsequently releases the director's discretionary early release science (ERS) data that include NIRCam, NIRSpec, and NIRISS data taken by the ERS programs of Cosmic Evolution Early Release Science (CEERS; Finkelstein et al. 2017, 2022c) and GLASS James Webb Space Telescope Early Release Science (GLASS; Treu et al. 2022). Further releases will deliver data of Cycle 1 observations that include Public Release IMaging for Extragalactic Research (PRIMER; Dunlop et al. 2021), UNCOVER (Labbe et al. 2021), and COSMOS-Webb (Kartaltepe et al. 2021) once the observations complete. Programs of guaranteed time observations (GTO), such as JWST Advanced Deep Extragalactic Survey (JADES; Bunker et al. 2020) will be also completed in the early years.  \nThis is a great development of observational astronomy, presenting the unprecedentedly deep and high quality data covering the infrared band ( > 2 µ m). In fact, after the releases of the ERO and ERS data sets, we find the explosive progresses of galaxy formation studies. The mass models of the ERO target cluster, SMACS J0723, are improved with the NIRCam imaging and NIRSpec spectroscopic data (Mahler et al. 2022; Pascale et al. 2022; Caminha et al. 2022). High redshift galaxies are searched in the ERO SMACS J0723 and ERS CEERS and GLASS fields, and are identified at z ∼ 9 -20 (Naidu et al. 2022b; Castellano et al. 2022; Leethochawalit et al. 2022; Adams et al. 2022; Morishita & Stiavelli 2022; Atek et al. 2022; Donnan et al. 2022; Finkelstein et al. 2022b; Yan et al. 2022). The morphological properties are investigated with the NIRCam images of the ERO SMACS J0723 and the CEERS observations via the comparisons of HST images for galaxies at z ∼ 3 -6 (Ferreira et al. 2022) and the rest-frame optical and near-infrared bands for galaxies at z ∼ 1 -2 (Suess et al. 2022), respectively. The infrared photometric properties of galaxies at z ∼ 1 -2 are studied with the NIRCam and MIRI images of the ERO SMACS J0723  \nobservations in conjunction with the ALMA archival data (Cheng et al. 2022). The ERO NIRSpec observations in SMACS J0723 provide high-quality spectra that allow to identify 10 galaxies at z = 1 . 2 -8 . 5, three of which reside at z = 7 . 7 -8 . 5 (Carnall et al. 2022), and to characterize the inter-stellar medium of the galaxies (Schaerer et al. 2022; Curti et al. 2022). NIRISS spectroscopic data complements the NIRSpec observations and provide spectroscopic sample of z ∼ 1 -8 galaxies (Roberts-Borsani et al. 2022; Wang et al. 2022; Boyett et al. 2022; Marchesini et al. 2022). More JWST results for galaxy formation are being actively reported.  \nIn this paper, we present a comprehensive study on high redshift galaxies using the first JWST/NIRCam datasets taken by the ERO and ERS programs. The deep infrared imaging data taken with NIRCam allow us to search for galaxies at z glyph[greaterorsimilar] 9, and to constrain the UV luminosity function and the cosmic SFR density in the universe 600 Myrs after the Big Bang. We will also perform thorough comparisons of galaxies found in our work and recent JWST studies.  \nThis paper is organized as follows. Section 2 presents the JWST observational data sets used in this study. Section 3 explains our sample selection and galaxy photometry catalog. In Section 4, we describe the mass model for the lensing cluster. We show our main results of UV luminosity functions and cosmic SFR densities in Section 5, and discuss the physical properties of early galaxies in Section 6. Section 7 summarizes our findings. Throughout this paper, we use the Planck cosmological parameter sets of the TT, TE, EE+lowP+lensing+BAO result (Planck Collaboration et al. 2020): Ω m = 0 . 3111, Ω Λ = 0 . 6899, Ω b = 0 . 0489, h = 0 . 6766, and σ 8 = 0 . 8102. All magnitudes are in the AB system (Oke & Gunn 1983).", '2.1. JWST/NIRCam Data': "We use four JWST NIRCam datasets obtained in the ERO and ERS programs, ERO SMACS J0723, ERO Stephan's Quintet, ERS CEERS, and ERS GLASS (Table 1). The total area is ∼ 90 arcmin 2 . We retrieved law data ( uncal.fits) from the MAST archive and reduced the data using the JWST pipeline version 1.6.3 development version (1.6.3.dev34+g6889f49, almost the same as 1.7.0). We use the Calibration Reference Data System (CRDS) context file of jwst 0995.pmap released in October, whose calibration values were derived using calibration observations of three different standard stars placed in all of the 10 NIRCam detectors. These new flux calibrations were verified using imaging of the globular cluster M92 (Boyer et al. 2022). In addition  \nTable 1. Limiting Magnitudes of the JWST Data  \nNote -Columns: (1) Field. (2) Effective area in arcmin 2 . (3)-(12) Typical limiting magnitudes which correspond to 5 σ variations in the sky flux measured with a circular aperture of 0 . '' 2-diameter in the deepest region.  \nFigure 1. Comparison of magnitudes. Magnitudes measured in the JWST F 150 W -band (left) and F 356 W -band (right) are compared with those in the HST F 160 W -band and Spitzer [3.6]-band, respectively. The measured magnitudes agree well with those in the HST and Spitzer images within ∼ 10%, indicating that the flux is reasonably calibrated. Note that we include a 10% error floor on all measured fluxes to account for possible systematic uncertainties.  \n<!-- image -->  \n21  \n22  \n23  \n24  \n25  \nF356W (JWST)  \nusing iraf tasks geomap and geotran . To check the reliability of the flux calibration, we compare our measured magnitudes in the JWST images with those in the HST and Spitzer images. As shown in Figure 1, the measured fluxes are almost consistent with those in the HST and Spitzer images, indicating that the flux is reasonably calibrated.  \nThe limiting magnitudes were measured in 0 . '' 1, 0 . '' 2, and 0 . '' 3-diameter circular apertures by randomly placing apertures in sky areas using Python packages Astropy/photutils . Sky areas were defined as pixels without objects detected by SExtractor . We measured the limiting magnitudes in bins of the weight values to take into account inhomogeneity of the depth. Us-  \n25  \n24  \n23  \n22  \n21  \nto the standard reduction, we added some processes to obtain better reduced images as follows. Before the Stage 2 calibration, we subtracted stray light features called 'wisps' by using a script provided by the NIRCam team 1 , and removed striping by using a script provided in the CEERS team (Bagley et al. 2022b) 2 . We ran the SkyMatch step individually on each frame of Stage 2 calibrated data before Stage 3 calibration, following a suggestion by the CEERS team (Bagley et al. 2022b). The images were pixel-aligned with a pixel scale of 0 . '' 015 / pixel, except for ones in the Stephan's Quintet field with a scale of 0 . '' 03 / pixel to reduce image sizes. Because the pipeline-processed images still showed a gradient of the sky background, we further subtracted the sky background using SExtractor (Bertin & Arnouts 1996). The intracluster light around the cluster center of SMACS J0723 is also removed in this process, although more sophisticated processes are sometimes employed (e.g., Livermore et al. 2017; Pagul et al. 2021). Note that our galaxy samples are not severely affected by systematics due to the intracluster light removal, because as shown in Section 3.5, the only candidate selected in the SMACS J0723 field, SM-z12-1, is located in the parallel field. The results of luminosity functions do not change beyond the errors if we remove the datapoint estimated with the SMACS J0723 data. Finally, we corrected for an astrometric offset between each detector and band  \n[3.6] (Spitzer)  \ning the weight map, we masked some regions around the edge of the detectors whose exposure time is short. We also measured a full-widths at half-maximum (FWHM) of the point spread function (PSF) in each image by selecting stellar objects in the magnitude-FWHM diagram. The measured limiting magnitudes in a 0 . '' 2diameter circular aperture, effective areas, and typical FWHMs of the PSFs are presented in Table 1. Here the effective area is defined as an area that is observed with all available bands before the foreground removal. The effect of the foreground will be taken into account in the completeness estimate (Section 5.1). In the following sections we detail our observational dataset in each field.", '2.1.1. ERO: SMACS J0723': 'A massive galaxy cluster at z = 0 . 39, SMACS J0723, was deeply observed with NIRCam, NIRSpec, MIRISS, and MIRI in the ERO (ERO-2736). The NIRCam images were taken in the six bands of F 090 W , F 150 W , F 200 W , F 277 W , F 356 W , and F 444 W , covering 11 . 0 arcmin 2 . The exposure time in each filter is ∼ 7500 seconds, and the 5 σ limiting magnitude in the F 356 W band is 29.9 mag.', "2.1.2. ERO: Stephan's Quintet": "Stephan's Quintet, a group of five local galaxies, was observed with NIRCam and MIRI in the ERO (ERO2732). The NIRCam images were taken in the six bands of F 090 W , F 150 W , F 200 W , F 277 W , F 356 W , and F 444 W , covering 42 arcmin 2 . We have masked central regions of the field that is affected by the five local galaxies in Stephan's Quintet, resulting in an effective area of 37 . 2 arcmin 2 , corresponding to ∼ 4 NIRCam pointings. The exposure time in each filter is roughly ∼ 1200 seconds, and the 5 σ limiting magnitude in the F 356 W band is 28.9 mag.", '2.1.3. ERS: CEERS': 'A part of the HST/CANDELS Extended Groth Strip (EGS) field is observed with JWST in the Cosmic Evolution Early Release Science (CEERS) survey (ERS-1345; Finkelstein et al. 2017, Finkelstein et al. 2022c). We use four pointing datasets of NIRCam obtained in June 2022 with the seven bands of F 115 W , F 150 W , F 200 W , F 277 W , F 356 W , F 410 M , and F 444 W , covering a total of 33 arcmin 2 . The exposure time in each filter is ∼ 2800 -6200 seconds, and the 5 σ limiting magnitude in the F 356 W band is 29.7 mag. Since the exposure times are not uniform across the four NIRCam pointings, we separately analyze the four pointing data.  \nA massive galaxy cluster, Abell 2744, was observed with JWST in the ERS program of Through the Looking GLASS (ERS-1324; Treu et al. 2017, 2022). Deep NIRCam images were taken in June 2022 in a parallel mode of NIRISS observations targeting the center of the cluster, in seven bands F 090 W , F 115 W , F 150 W , F 200 W , F 277 W , F 356 W , and F 444 W , covering 6 . 8 arcmin 2 . The exposure time in each filter is ∼ 5600 -23000 seconds, and the 5 σ limiting magnitude in the F 356 W band is 29.9 mag. The lensing magnification is negligible in the field of the NIRCam observations that is ∼ 5 arcmin away from the cluster center of Abell 2744.', '2.2. JWST/NIRSpec': 'We use publicly available data from the ERO NIRSpec observations targeting the field of the SMACS J0723 cluster (ERO-2736). The NIRSpec observations consist of two pointings with the same Multi Shutter Array (MSA) configuration. NIRSpec observations were carried out by using the disperser-filter combinations of G235M/F170LP and G395M/F290LP, which cover the wavelength range from 1.66 µ m to 5.16 µ m with a spectral resolution of R ∼ 1000. The total exposure time of the two individual pointings is 8,840 seconds for each grating. The NIRSpec observations have taken spectra for a total of 35 objects, three of which, s04590 ( z spec = 8 . 495), s06355 ( z spec = 7 . 664), and s10612 ( z spec = 7 . 659), are securely identified at z > 7 whose dropouts can be covered with the NIRCam F 090 W band. See K. Nakajima et al. in prep. for details of the reduction and analysis of the NIRSpec data.', '2.3. HST/ACS and WFC3': "HST multi-band images are available in the fields of SMACS J0723 and CEERS (EGS). We downloaded HST ACS and WFC3 images in the SMACS J0723 and CEERS fields from the websites of RELICS (Coe et al. 2019) 3 and CEERS 4 , respectively. We found a small offset ( ∼ 0 . '' 2) of the WCS between the JWST and HST images. In this paper we use coordinates of the JWST images.", '3.1. Photometric Catalog': "We construct multi-band source catalogs from the JWST data to select the F 115 W , F 150 W , and F 200 W -dropout galaxies. We use SWarp (Bertin et al. 2002)  \nto produce our detection image that is a weighted mean image of the bands redder than the Lyman break in each dropout selection (i.e., F 150 W , F 200 W , F 277 W , F 356 W , F 410 M , and F 444 W for the F 115 W -dropout selection, F 200 W , F 277 W , F 356 W , F 410 M , and F 444 W for the F 150 W -dropout selection, and F 277 W , F 356 W , F 410 M , and F 444 W for the F 200 W -dropout selection). To measure object colors, we match the image PSFs to the F444W-band images whose typical FWHM of the PSF is glyph[similarequal] 0 . '' 16, the largest of the NIRCam multi-band images.  \nWe perform source detection and photometry with SExtractor (version 2.5.0; Bertin & Arnouts 1996). We found that the photometry with SExtactor MAG AUTO performs better as a total magnitude than that with photutils isophotal flux , which is one of the outputs of the JWST calibration pipeline, using CEERS simulated images with a mock galaxy catalog created with the Santa Cruz Semi-Analytic Model (Somerville et al. 2021; Yung et al. 2022). 5 We run SExtractor in dual-image mode for each image with its detection image, having the parameter set as follows: DETECT MINAREA = 5, DETECT THRESH = 3 . 0, ANALYSIS THRESH = 3 . 0, DEBLEND NTHRESH = 32, and DEBLEND MINCOUNT = 0 . 005. The total number of the objects detected is ∼ 250,000. We measure the object colors with the MAG APER magnitudes defined in a 0 . '' 3diameter circular aperture in the PSF-matched images. Source detections are evaluated with 0 . '' 1 and/or 0 . '' 2diameter circular apertures in the original (not PSFmatched) images. The total magnitudes are estimated from the 0 . '' 3-diameter aperture magnitudes with the aperture correction. The value of the aperture correction is defined as the difference between the MAG AUTO magnitude and the 0 . '' 3-diameter aperture magnitude in an image of the weighted mean of the PSF-matched images whose wavelengths are longer than the Lyman break (i.e., F 150 W , F 200 W , F 277 W , F 356 W , F 410 M , and F 444 W for the F 115 W -dropout selection, F 200 W , F 277 W , F 356 W , F 410 M , and F 444 W for the F 150 W -dropout selection, and F 277 W , F 356 W , F 410 M , and F 444 W for the F 200 W -dropout selection). Furthermore, we correct for a small offset ( ∼ 0 . 1 mag) between the measurement of MAG AUTO and the true total magnitude due to the wing of the PSF not captured with MAG AUTO (see Sections 2.2 and 2.5.1 in Finkelstein et al. 2022a). We measure this offset by inserting mock galaxy images randomly in the real images, and measure the magnitudes using SExtractor in a similar manner  \nas the completeness simulation described later in Section 5.1. To account for systematic uncertainties of the flux measurements (e.g., zero-point correction), we include 10% error floor on all measured fluxes. Finally we correct for the galactic extinction using Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) and make final photometric catalogs.", '3.2. Dropout Selection': "From the photometric catalogs constructed in Section 3.1, we construct z ∼ 9 -16 dropout galaxy catalogs based on the Lyman break color selection technique (e.g., Steidel et al. 1996; Giavalisco 2002). As shown in Figure 2, galaxy candidates at z ∼ 9, 12, and 16 can be selected by the F 115 W , F 150 W , and F 200 W -dropout selections, respectively.  \nFirst, to identify secure sources, we select source whose signal-to-noise ratios (S/Ns) in a 0 . '' 2-diameter circular aperture are higher than 5 in the detection images. We also require sources to be detected at > 3 . 5 σ levels in at least two bands redder than the Lyman break. We then select dropout galaxy candidates by using their broadband spectral energy distribution (SED) colors. We adopt the following color criteria:  \nF 115 W -dropout ( z ∼ 9):  \n( F 115 W -F 150 W > 1 . 0) ∧ (1)  \n( F 150 W -F 277 W < 1 . 0) ∧ (2)  \n( F 115 W -F 150 W > ( F 150 W -F 277 W ) + 1 . 0) (3)  \nF 150 W -dropout ( z ∼ 12):  \n( F 150 W -F 200 W > 1 . 0) ∧ (4)  \n( F 200 W -F 356 W < 1 . 0) ∧ (5)  \n( F 150 W -F 200 W > ( F 200 W -F 356 W ) + 1 . 0) (6)  \nF 200 W -dropout ( z ∼ 16):  \n( F 200 W -F 277 W > 1 . 0) ∧ (7)  \n( F 277 W -F 444 W < 1 . 0) ∧ (8)  \n( F 200 W -F 277 W > 1 . 5( F 277 W -F 444 W ) + 1 . 0) (9)  \nWe select sources with prominent breaks with the criteria of Equations (1), (4), and (7), and measure the slope of the continuum and remove red interlopers with Equations of (2-3), (5-6), and (8-9). To measure the slope of the continuum, we use the bands that have a large wavelength difference as much as possible, while not biting the Balmer break in the redder band. To remove foreground interlopers, we exclude  \nF150W -F277W --Figure 2. Two-color diagrams of F 115 W -F 150 W vs. F 150 W -F 277 W (left), F 150 W -F 200 W vs. F 200 W -F 356 W (center), and F 200 W -F 277 W vs. F 277 W -F 444 W (right) corresponding to the color selections for F 115 W -dropouts at z ∼ 9, F 150 W -dropouts at z ∼ 12, and F 200 W -dropouts at z ∼ 16, respectively. The red squares represent our dropout galaxy candidates that meet the color selection criteria indicated with the red lines. The blue lines denote colors of the dropout galaxy models with UV spectral slopes of β UV = -2 . 3 and -1 . 3 whose redshifts are indicated with the numbers and the blue circles with an interval of ∆ z = 0 . 2. The black dotted, dashed, dot-dashed lines show colors of typical elliptical, Sbc, and irregular galaxies (Coleman et al. 1980) redshifted from z = 0 to 7. The star marks present expected colors of Galactic dwarf stars (Patten et al. 2006; Kirkpatrick et al. 2011). These colors of dwarf stars are estimated by interpolating available flux measurements obtained by ground based telescopes ( J , H , and K -bands) and the Spitzer telescope ([3.6] and [4.5]).  \n<!-- image -->  \nsources with continuum detections at > 2 σ levels in the 0 . '' 1- or 0 . '' 2-diameter apertures in bands bluer than the Lyman break, i.e., the F 090 W band for the F 115 W -dropouts, F 090 W and F 115 W bands for the F 150 W -dropouts, and F 090 W , F 115 W , and F 150 W bands for the F200W-dropouts. To select reliable candidates, we restrict our dropout selections in fields where the bluer band than the Lyman break is available; the F 115 W -dropout selection is only performed in the GLASS field. We also apply a criterion of SExtractor stellarity parameter, CLASS STAR , of < 0 . 9, to remove stellar contaminants. Finally, we visually inspect images of the selected sources to remove spurious sources or sources affected by nearby bright objects and diffraction spikes of bright stars. We removed about half of the selected objects in this process. We also visually inspect HST images of the selected sources in the SMACS J0723 and CEERS fields to check whether the source is consistent with being a high redshift galaxy, although the HST images are typically ∼ 1 -2 mag shallower than the JWST images in these fields.", '3.3. SED Fitting': 'To further remove low redshift interlopers, we perform galaxy SED fitting with the flexible Bayesian inference code prospector (Johnson et al. 2021), and derive the photometric redshift. Model spectra are derived from Flexible Stellar Population Synthesis (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) package with the modules for Experiments in Stellar Astrophysics Isochrones  \nand Stellar Tracks (MIST; Choi et al. 2016). The boost of ionizing flux production of massive stars are included in the MIST isochrones (Choi et al. 2017). Here we assume the stellar initial mass function (IMF) determined by Chabrier (2003), the Calzetti et al. (2000) dust extinction law, and the intergalactic medium (IGM) attenuation model by Madau (1995). Note that the choice of the IGM attenuation model does not affect our galaxy selection at z ∼ 9 -16 because the flux bluer than the Ly α break is almost absorbed by the highly neutral IGM at these redshifts regardless of the choice of the IGM attenuation model. The Ly α emission line is also masked considering the high IGM neutral fraction at these redshifts. We adopt a flexible star formation history with five bins that are spaced equally in logarithmic times between 0 Myr and a lookback time that corresponds to z = 30, where the SFR within each bin is constant. We change the redshift, optical depth in the V -band, star formation history, and total stellar mass as free parameters, while fix the metallicity to Z = 0 . 2 Z glyph[circledot] . We assume a continuity prior for the star formation history, and flat priors for other parameters in the range 0 < z < 20, 0 < τ V < 2, and 6 < log( M ∗ /M glyph[circledot] ) < 12. We search for the best-fit model to the observed photometry with the MCMC method by using emcee (Foreman-Mackey et al. 2013).  \nBased on the results of the SED fitting, we select objects whose high-redshift solution is more likely than the low-redshift ones, by using ∆ χ 2 , which is defined as a difference between the χ 2 values of the best high-  \nhift solution and the lower-redshift solution, ∆ χ 2 = χ 2 ( z low ) -χ 2 ( z high ). Previous studies use a criterion of ∆ χ 2 > 4 . 0, corresponding to a 2 σ level (e.g., Bowler et al. 2020; Harikane et al. 2022a; Donnan et al. 2022; Finkelstein et al. 2022b). However, given the small number of the available JWST bands bluer than the Lyman break and the expected small number density of z > 9 galaxies, it is possible that this criterion is not sufficient to remove low redshift interlopers. To determine the threshold value for ∆ χ 2 , we use the CEERS simulated NIRCam images. In the CEERS simulated images, mock galaxies at z = 0 -10 in Yung et al. (2019, 2022) are inserted using the JWST data simulator Mirage (Hilbert et al. 2019). We measure fluxes of mock galaxies in each band in the same manner as our real dropout galaxy selection (Section 3.1), and conduct the SED fitting using prospector . As shown in Figure 3, at least eight sources at z true ∼ 3 -4 in the simulations have the best photometric redshifts of z phot ∼ 12 -15 and ∆ χ 2 = 4 -9, indicating that the criterion of ∆ χ 2 > 4 is not sufficient to remove low redshift interlopers. Thus in this study, we instead adopt a strict screening criterion of ∆ χ 2 > 9 . 0, which can remove these interlopers. The inclusion of this strict criterion does not introduce a bias with respect to a color of the UV continuum for bright galaxies, because the strength of the break is the most important factor to determine the ∆ χ 2 value. For faint galaxies, about 40% of them at ∼ 29 -30 mag will be missed due to the inclusion of this criterion, and this effect is taken into account in the completeness estimate in Section 5.1.', '3.4. Comparisons with Spectroscopic Redshifts': 'To test the reliability of our galaxy selections and SED fitting, we compare our photometric redshift estimates with the spectroscopic results obtained in the NIRSpec observations (Section 2.2). Since there are currently no z > 9 source spectroscopically confirmed with NIRSpec in the filelds used in this study, we focus on the three galaxies at z > 7, s04590 ( z spec = 8 . 495), s06355 ( z spec = 7 . 664), and s10612 ( z spec = 7 . 659). We measure fluxes of the three galaxies in the NIRCam images and estimate photometric redshifts using prospector , in the same manner as our dropout galaxies. Figures 4 and 5 present results of the SED fitting and comparison with the spectroscopic redshifts. We found that the estimated photometric redshifts agree well with the spectroscopic redshifts within ∼ 2 σ uncertainties, indicating that our SED fitting works well to estimate the redshift from the NIRCam photometry.', '3.5. Final Sample': "Finally we select 13, 8, and 2 dropout galaxy candidates at z ∼ 9, 12, and 16, respectively (Table 2). The photometric redshifts range from z ∼ 8 . 7 to 16 . 4, demonstrating the power of JWST exploring the early universe (Figure 6). Examples of the snapshots and SEDs of selected galaxies at z ∼ 9, 12, and 16 are presented in Figures 7 and 8. These sources show a sharp discontinuity around the Lyman break band, a flat or blue continuum, and non-detection in the bluer bands than the break, all of which are consistent with a high redshift galaxy. Photometric properties of our galaxy candidates are summarizes in Tables 3-5. Note that no objects appear in more than one final dropout sample.  \nTo investigate morphological properties of our galaxy candidates, we fit our galaxy candidates with the S'ersic profile using galfit (Peng et al. 2010). We find that some bright candidates, i.e., GL-z9-1, GL-z12-1, CR2z12-1, CR-z16-1, and S5-z16-1, are clearly more extended than the PSF, although GL-z12-1 is compact compared to other candidates, implying potential AGN activity. We stack images of other faint candidates at each redshift, and find that the stacked images also show extended profiles with respect to the PSF. We thus conclude that the stellar contamination is negligible. Details of the morphological properties of our candidates are presented in Ono et al. (2022).  \nOne of the highest redshift source candidates in our catalogs is CR2-z16-1 at z = 16 . 25 +0 . 24 -0 . 46 in the CEERS2 field. As discussed later in Section 3.6.6, CR2-z16-1 is firstly identified as a z = 16 . 4 source (ID 93116) in Donnan et al. (2022). As shown in the middle panel of Figure 8, our measured fluxes are almost consistent with those presented in Naidu et al. (2022a) and Finkelstein et al. (2022c), while fluxes in Donnan et al. (2022) are systematically fainter than our measurements, probably because Donnan et al. (2022) assume the PSF for the aperture correction. The colors measured in these three studies (this study, Donnan et al. 2022, Naidu et al. 2022a, and Finkelstein et al. 2022c) consistently show a clear break around F 200 W -band, consistent with a z = 16 . 3 galaxy, although there are some discussions about a possible solution of a dusty line emitter at z ∼ 5 (Zavala et al. 2022; Naidu et al. 2022a). Although the NIRSpec spectroscopy is required to conclude the redshift of CR2-z16-1, we include this source as a F 200 W -dropout galaxy candidate.  \nThe other candidate at z ∼ 17 is S5-z16-1 at z = 16 . 41 +0 . 66 -0 . 55 identified in the Stephan's quintet field. Although this source is located in a region whose exposure time is relatively short compared to the central region of the field, our position-dependent estimates of flux uncertainties indicate that the source detection, color, and", 'Interlopers in Simulations': "Figure 3. Examples of mock galaxies whose true redshifts are z true ∼ 3 but selected as F 150 W -dropout galaxies at z phot ∼ 12 -15 identified in the simulated NIRCam images. For each object, the top left panel shows the 1 . '' 5 × 1 . '' 5 snapshots in NIRCam bands with a 3-pixel smoothing whose band names are indicated with the red labels. The bottom left panel presents the SED of the object. The red symbols with error bars are measured magnitudes or 2 σ upper limits, and the blue curve shows the best-fit model. The true and estimated redshifts with 2 σ errors are indicated with the black and blue texts, respectively. The χ 2 value is shown in the right panel as a function of redshift. The black curve is the true SEDs at z ∼ 3, while the blue curve denotes χ 2 values of our SED fitting for our photometric redshift determination. These objects meet the weak photometric redshift criterion of ∆ χ 2 > 4, but do not meet our strict criterion of ∆ χ 2 > 9, where ∆ χ 2 is the χ 2 difference between the best high redshift solution and a lower redshift solution, ∆ χ 2 = χ 2 ( z low ) -χ 2 ( z high ). See texts for details.  \n<!-- image -->  \nFigure 4. Top: NIRSpec spectrum of s04590 at a spectroscopic redshift of z = 8 . 495. The spectroscopic redshift is confirmed with the H β , H γ , H δ , [Oiii] λλ 5007,4959, [Oiii] λ 4363, [Oii] λ 3727, [Ne iii ] λ 3967, [Ne iii ] λ 3869, and the tentative Ciii] λ 1909 lines. The flux is arbitrary. Bottom left: Optical to near-infrared SED of s04590. The red circles and arrows indicate the measured magnitudes and 2 σ limits, respectively. The filled (open) red symbols denote the measurements and the limits obtained with JWST/NIRCam (HST/ACS and WFC3). The blue curve and the blue redshift label represent the best-fit model SED and the photometric redshift with 2 σ errors derived by our photometric redshift technique, which is compared with the spectroscopic redshift indicated with the black label. The images on this panel show 1 . '' 5 × 1 . '' 5 cutout images of s04590 in the NIRCam bands with a 3-pixel smoothing whose band names are indicated with the red labels. Bottom right: χ 2 as a function of redshift. The blue curve denotes χ 2 values of our SED fitting for our photometric redshift determination. The vertical dashed line indicates the spectroscopic redshift that agrees well with the photometric redshift.  \n<!-- image -->  \nFigure 5. Same as Figure 4, but for s06355 at z spec = 7 . 664 and s10612 at z spec = 7 . 659. Our estimates of the photometric redshifts agree well with the spectroscopic redshifts.  \n<!-- image -->  \nFigure 6. Absolute UV magnitude as a function of redshift for galaxies at 6 < z < 16. The red diamonds represent our dropout galaxy candidates selected with the JWST images. The red open circles show HD1 and HD2 previously found by the combination of the images taken with Sptizer and ground-based telescopes (Harikane et al. 2022a). The gray square and circles denote GN-z11 (Oesch et al. 2016; Jiang et al. 2021) and dropout galaxies selected with deep HST images (Bouwens et al. 2015).  \n<!-- image -->  \nTable 2. Number of Our Dropout Candidates  \nnon-detections are robust against the uncertainties. An emission line feature is detected with ALMA in S5-z161 (Fujimoto et al. 2022), which would be interpreted as either [Oiii] 52 µ m at z = 16 . 01 or [Cii]1 58 µ m at z = 4 . 61. We include this possible candidate in the luminosity function calculation, although the luminosity is remarkably high compared to our expectations at this high redshift.  \nFigure 7. (Top:) The left panel presents the optical to near-infrared SEDs of the z ∼ 9 dropout galaxy, GL-z9-1. The red circles and arrows show the measured magnitudes and 2 σ upper limits, respectively. The blue curve denotes the best-fit model SED whose redshift and 2 σ errors are presented in the upper left with the blue labels. The gray curve is a significantly worse fit of a low redshift solution. The images on this panel are 1 . '' 5 × 1 . '' 5 cutout images of GL-z9-1 in the NIRCam bands with a 3-pixel smoothing whose band names are indicated with the red labels. The right panel shows χ 2 values of the SED fitting as a function of redshift. (Middle:) Same as the top panels, but for another z ∼ 9 dropout galaxy candidate, GL-z9-2. (Bottom:) Same as the top panels, but for a z ∼ 12 dropout galaxy candidate, GL-z12-1.  \n<!-- image -->  \nFigure 8. Same as Figure 7, but for z ∼ 12 dropout galaxy candidate, CR2-z12-1 (top), and z ∼ 16 dropout galaxy candidates, CR2-z16-1 (middle) and S5-z16-1 (bottom). The open red symbols denote the measurements and the limits obtained with HST/ACS and WFC3. The orange, pink, and brown open symbols in the middle panel are measurements in Donnan et al. (2022), Naidu et al. (2022a), and Finkelstein et al. (2022c), respectively.  \n<!-- image -->  \nTable 3. Summary of Our F 115 W -Dropoout Candidates at z ∼ 9  \nNote -(1) Name. (2) Right ascension. (3) Declination. (4) Total magnitude in the F 356 W band with 1 σ errors. (5) F 115 W -F 150 W color with 1 σ errors. (6) F 150 W -F 277 W color with 1 σ errors. (7) Absolute UV magnitude with 1 σ errors. (8) Photometric redshift with 2 σ errors. (9) χ 2 difference between the best high redshift solution and a lower redshift solution, ∆ χ 2 = χ 2 ( z low ) -χ 2 ( z high ). (10) Reference (This: this work, N22: Naidu et al. (2022b), C22: Castellano et al. (2022), D22: Donnan et al. (2022)) and note for a reason why the source is not selected in this study (1: > 2 σ detection in F 090 W , 2: ∆ χ 2 < 9 . 0).  \nTable 4 . Summary of Our F 150 W -Dropoout Candidates at z ∼ 12  \nTable 4 continued  \nTable 4 (continued)  \nTable 4 continued  \nJWST Galaxies at z ∼ 9 -17  \nTable 4 (continued)  \nNote -(1) Name. (2) Right ascension. (3) Declination. (4) Total magnitude in the F 356 W band with 1 σ errors. (5) F 150 W -F 200 W color with 1 σ errors. (6) F 200 W -F 356 W color with 1 σ errors. (7) Absolute UV magnitude with 1 σ errors. Values of galaxies in the SMACS J0723 field are after the lensing magnification correction with glafic . (8) Photometric redshift with 2 σ errors. (9) χ 2 difference between the best high redshift solution and a lower redshift solution, ∆ χ 2 = χ 2 ( z low ) -χ 2 ( z high ). (10) Reference (This: this work, N22: Naidu et al. (2022b), C22: Castellano et al. (2022), D22: Donnan et al. (2022), F22: Finkelstein et al. (2022b), Y22: Yan et al. (2022)) and note for a reason not selected in this study (1: > 2 σ detection in F 090 W , 2: F 150 W -F 200 W < 1 . 0, 3: ∆ χ 2 < 9 . 0, 4: < 5 σ in the detection image).  \nTable 5 . Summary of Our F 200 W -Dropoout Candidates at z ∼ 16  \nNote -(1) Name. (2) Right ascension. (3) Declination. (4) Total magnitude in the F 356 W band with 1 σ errors. (5) F 200 W -F 277 W color with 1 σ errors. (6) F 277 W -F 444 W color with 1 σ errors. (7) Absolute UV magnitude with 1 σ errors. Values of galaxies in the SMACS J0723 field are after the lensing magnification correction with glafic . (8) Photometric redshift with 2 σ errors. (9) χ 2 difference between the best high redshift solution and a lower redshift solution, ∆ χ 2 = χ 2 ( z low ) -χ 2 ( z high ). (10) Reference (This: this work, D22: Donnan et al. (2022)) and note for a reason not selected in this study (1: F 200 W -F 277 W < 1 . 0, 2: ∆ χ 2 < 9 . 0).", '3.6. Comparison with Previous Studies': 'Some other studies identified galaxy candidates at z > 9 using the JWST NIRCam ERO and/or ERS  \ndatasets. Here we review these studies and compare their samples with our galaxy samples. Tables 3-5 summarize properties of other possible candidates that were  \nselected in other studies but did not meet our selection criteria. These comparisons were conducted on November 20, 2022, and we clarify the version of the paper we compared in the following sections.', '3.6.1. Naidu et al. (2022b)': 'Using the ERS CEERS and GLASS datasets, Naidu et al. (2022b) found two bright galaxy candidates at z ∼ 10 and 12, GLASS-z10 and GLASS-z12, which correspond to GL-z9-1 and GL-z12-1 in our sample, respectively. Their estimates of the photometric redshifts ( z = 10 . 35 +0 . 38 -0 . 51 and z = 12 . 38 +0 . 13 -0 . 27 for GL-z101 and GL-z12-1, respectively, with prospector , from the ApJL published version) are consistent with our estimates ( z = 10 . 49 +0 . 53 -0 . 72 and z = 12 . 28 +0 . 08 -0 . 07 ).', '3.6.2. Castellano et al. (2022)': 'Castellano et al. (2022) identified seven galaxy candidates at z ∼ 9 -12 with the color selection using the ERS GLASS dataset. Among the six candidates from the ApJL published version, three candidates, GHZ1, GHZ2, and GHZ4, are selected in our selection. GHZ1 (GHZ2) is GL-z9-1 (GL-z13-1) in our sample, and their photometric redshift, z = 10 . 53 -10 . 63 ( z = 12 . 11 -12 . 30) is comparable with our estimates. GHZ4 was identified in our selection as GL-z9-2, and their photometric redshift ( z = 9 . 93 -10 . 08) agrees with our estimate ( z = 10 . 46 +0 . 45 -0 . 99 ). The other three candidates, GHZ3, GHZ5, and GHZ6, did not meet our selection criteria, due to a possible detection in the F 090 W -band or ∆ χ 2 < 9, although GHZ3 and GHZ5 were also reported to have low redshift solutions in Castellano et al. (2022).', '3.6.3. Leethochawalit et al. (2022)': 'Leethochawalit et al. (2022) studied photometric properties of galaxies at 7 < z < 9 using the ERS GLASS dataset. We refer to the manuscript version 2 that was submitted to arXiv on October 4. Since their galaxies are identified from the F 090 W -dropout selection, we do not expect significant overlap between their galaxy sample and ours.', '3.6.4. Adams et al. (2022)': 'Adams et al. (2022) identified four galaxy candidates at 9 < z < 12 in the ERO dataset taken in the SMACS J0723 field. We refer to the manuscript version 2 that was submitted to arXiv on August 9. Since three of them are expected to be at z < 10, it is reasonable that they are not identified in our F 150 W -dropout selection ( z ∼ 12). The other source (ID 10234) is estimated to be at z = 11 . 42, around the edge of our redshift selection window (see Figure 11), and not selected in  \nour study due to its insufficient F 150 W -F 200 W color and a possible detection in the F 090 W -band.', '3.6.5. Atek et al. (2022)': 'Using the ERO dataset in the SMACS J0723 field, Atek et al. (2022) selected 10 galaxy candidates at 10 < z < 16. We refer to the manuscript version 2 that was submitted to arXiv on October 31. Among them, four candidates, SMACS z12a, SMAC z12b, SMACS z16a, and SMACS z16b, have photometric redshifts of z > 12, and are expected to be overlapped in our galaxy catalogs. However, none of them are selected as high redshift galaxy candidates in our study, due to their insufficient colors or ∆ χ 2 < 9.', '3.6.6. Donnan et al. (2022)': 'Donnan et al. (2022) selected 45 galaxies at z > 8 . 5 using ERO SMACS J0723 and ERS GLASS and CEERS datasets. We refer to the manuscript version 2 that was submitted to arXiv on October 22. Among the 45 galaxies, three galaxies (IDs 1698, 6415, and 17487) are identified in the GLASS dataset, and are also selected in this study as GL-z9-1, GL-z9-4, and GL-z12-1, respectively. Their photometric redshifts ( z = 10 . 45 +0 . 26 -0 . 16 , z = 10 . 79 +0 . 45 -0 . 66 , and z = 12 . 42 +0 . 27 -0 . 21 for IDs 1698, 6415, and 17487, respectively) are consistent with our estimates ( z = 10 . 49 +0 . 53 -0 . 72 , z = 10 . 19 +0 . 63 -0 . 55 , and z = 12 . 28 +0 . 08 -0 . 07 ). The brightest candidate in Donnan et al. (2022) is ID 93316 at z = 16 . 39 +0 . 32 -0 . 22 , which is CR2z16-1 at z = 16 . 25 +0 . 24 -0 . 46 in our catalog. In the version 2, Donnan et al. (2022) newly selected ID 32395 2 at z = 12 . 29 +0 . 91 -0 . 32 , which is also selected in this study as CR2-z12-1 at z = 11 . 63 +0 . 51 -0 . 53 , which was firstly identified in Finkelstein et al. (2022b) Donnan et al. (2022) presented other three candidates at z > 12, but these candidates are not selected in this study due to ∆ χ 2 < 9.', '3.6.7. Finkelstein et al. (2022b)': "One galaxy candidate at z ∼ 14, dubbed Maisie's Galaxy in Finkelstein et al. (2022b), is also selected in this study as CR2-z12-1. We refer to the manuscript version 2 that was submitted to arXiv on September 7. The photometric redshift presented in Finkelstein et al. (2022b) is z = 11 . 8 +0 . 3 -1 . 2 , consistent with our estimate ( z = 11 . 63 +0 . 51 -0 . 53 ).", '3.6.8. Yan et al. (2022)': 'Yan et al. (2022) identified a total of 88 galaxy candidates at z ∼ 11 -20 in the ERO SMACS J0723 field, including 63 F 150 W -dropouts and 15 F 200 W -dropouts, possibly overlapping with our F 150 W and F 200 W -dropout candidates, respectively. We refer to the manuscript version 1 that was submitted to arXiv  \nW  \n0  \n5  \n1  \nF  \n-  \nW  \n5  \n1  \n1  \nF  \n1  \n0  \n3  \n2  \n1  \n0  \n1  \n2  \n1  \n0  \n1  \n2  \nF277W  \nF444W  \nF150W -F277W --Figure 9. Same as Figure 2, but for evaluating the interlopers of foreground galaxies. The gray curves indicate colors of model galaxies at z = 0 -8 that are produced with PANHIT (Mawatari et al. 2020a). See texts for details of the models. The black arrow indicates a shift of the colors with dust extinction of ∆ E ( B -V ) = +0 . 1. The magenta circle in the right panel is a dusty starburst galaxy at z ∼ 5 that may appear as a z > 15 galaxy discussed in Zavala et al. (2022). Our color selection criteria avoid these low-redshift interlopers at z ∼ 0 -8.  \n<!-- image -->  \non July 23. Out of 61 and 15 sources in their F 150 W -dropouts and F 200 W -dropouts, we identify 54 and 11 objects in our original photometric catalogs, respectively. However, we cannot identify counterparts of the remaining 11 sources, F150DA-013, F150DA-047, F150DA-057, F150DB-004, F150DB-023, F150DB-056, F150DB-058, F200DB-015, F200DB-109, F200DB-175, and F200DB-181, probably because their SNRs are not sufficient to be identified in this study, the source is severely affected by nearby bright objects, or a WCS offset between Yan et al. (2022) and this study is too large to identify the counterparts. Among the 54 objects identified as F 150 W -dropouts, F150DA-053 at z = 11 . 71 +1 . 56 -0 . 54 is SM-z12-1 at z = 12 . 47 +1 . 19 -0 . 72 in this study. We have checked photometry of the other 53 and 11 objects identified in our original photometric catalogs, but none of them are selected as high redshift candidates in this study, due to their insufficient colors of the break, ∆ χ 2 < 9, and/or an insufficient S/N in the detection image.', '3.6.9. Summary of the Comparisons': 'In summary, we have found that bright candidates reported in previous studies are reproduced in this study, such as GL-z9-1, GL-z12-1, CR2-z12-1, and CR2-z161. However, some of faint candidates reported in other studies are not selected in our selection criteria, because most of these faint candidates are selected by photometric redshifts but with a weak criterion (e.g., ∆ χ 2 > 4) or by relatively weak color selection criteria (e.g., F 150 W -F 200 W > 0 . 5). It is expected that the contamination fraction in such faint candidates is high, given the small ∆ χ 2 values (see discussions in Section 3.3). These comparisons indicate that our selection  \ncriteria are conservative enough to remove foreground interlopers while keeping bright and reliable candidates.', '3.7. Contamination': 'We check whether our sample is largely contaminated by foreground interlopers or not. One of the major sources for contamination is low redshift galaxies whose Balmer breaks are redshifted to the wavelength of the Lyman break of our dropout galaxies. To test the effect of such contamination, we make a mock catalog of galaxies with Balmer breaks at z = 0 -8. We first generate model spectra of galaxies by using PANHIT (Mawatari et al. 2020a) assuming a delayedτ star formation history, τ = 0 . 01, 0 . 03, 0 . 1, 0 . 3, and 1 Gyr; stellar age of 0 . 01 -1 . 3 Gyr; and metallicity of Z = 0 . 0001, 0.0004, 0.004, 0.008, 0.02, and 0.05. In Figure 9 we plot model tracks of the z = 0 -8 galaxies with colors of our selected galaxy candidates and our color selection criteria. We find that the z = 0 -8 galaxies with Balmer breaks have a relatively small break color (e.g., F 150 W -F 200 W < 1 . 0 in the F 150 W -dropout selection) or larger break color and red continuum color ( F 150 W -F 200 W > 1 . 0 and F 200 W -F 356 W > 0). Our color selection criteria avoid model tracks of these z = 0 -8 galaxies, and most of our candidates are located far from these model tracks. In the right panel of Figure 9, we also plot a dusty starburst galaxy at z ∼ 5 that may appear as a z > 15 galaxy discussed in Zavala et al. (2022). Such a dusty interloper is also removed from our sample due to the red continuum color ( F 277 W -F 444 W > 1 . 0).  \nTo evaluate the effect of contamination from low redshift objects scattering into our selection criteria due  \nF115W-dropout  \n(  \nz  \n∼  \n9)  \nW  \n7  \n7  \n2  \nF  \n-  \nW  \n0  \n0  \n2  \nF  \n1  \n0  \n3  \n2  \nF200W-dropout  \n(  \nz  \n∼  \n16)  \nZavala+22', 'F115W-Dropout (stacked)': "F150W-Dropout (stacked)  \nFigure 10. Stacked images of our F 115 W -dropouts (top) and F 150 W -dropouts (bottom). The size of the images is 1 . '' 5 × 1 . '' 5. There are no positive signals found at the positions of the dropouts in the F 090 W image (top) and the F 090 W and F 115 W images (bottom) whose wavelength ranges (rest-frame < 1216 ˚ A) do not include emission from z ∼ 9 and 12 sources, indicating that our samples are not significantly contaminated by foreground interlopers.  \nto the photometric noise at the depth of the observations, we conduct Monte Carlo simulations using the real datasets in the same manner as previous studies (e.g., Bouwens et al. 2015; Ono et al. 2018; Harikane et al. 2022b). We start from multi-band catalogs constructed in the GLASS field whose images are sufficiently deep. We create 100 mock catalogs by perturbing the measured fluxes by adding photometric scatters based on the flux uncertainties in each band in the CEERS and Stephan's Quintet fields whose depths are shallower than the GLASS field. We select high redshift galaxies from the mock catalogs with the same selection criteria as our real selection. In the same manner as Bouwens et al. (2015), we classify sources that are selected but that show detections in the band bluer than the break in the original catalogs as contaminants. Based on these simulations, we find that the contamination rate due to the scatter is < 6% for the F 150 W and F 200 W -dropout selections. The contamination rates for the F 115 W -dropout selection cannot be evaluated with this procedure, because the galaxy selection is conducted only in the GLASS field where we start the simulation. However, the good agreement in z ∼ 9 number densities between our results and previous studies (Section 5.2) indicates that the contamination is not significant in our F 115 W -dropout sample.  \nTo further test the contamination, we stack images of the 13 F 115 W -dropout candidates and 8 F 150 W -dropout candidates. If the sample is significantly contaminated by low redshift interlopers, the stacked images should show signals in a band whose wavelength is bluer than the Lyman break. Figure 10 presents stacked images of our F 115 W - and F 150 W -dropout candidates. There are no significant positive signals found in the F 090 W -band and in the F 090 W - and F 115 W -bands for  \nthe F 115 W - and F 150 W -dropout candidates, respectively, suggesting that our samples are not significantly contaminated by low redshift interlopers. These tests and comparisons in Section 3.6 indicate that our conservative selection criteria with careful screening of low redshift interlopers provide a reliable sample of z ∼ 9 -16 galaxy candidates, suitable for statistical studies such as luminosity function measurements.", '4. MASS MODEL': "Various mass models for SMACS J0723 are produced by parametric mass modeling algorithms. The RELICS survey (Coe et al. 2019) team provides the mass models 6 developed with the Lenstool (Jullo et al. 2007; Fox et al. 2022) and glafic (Oguri 2010) codes, both of which are constructed using the HST data. A new model with the Lenstool code is constructed with the JWST ERO data in Mahler et al. (2022). Golubchik et al. (2022) recently present a mass model of SMACS J0723 developed by the light-traces-mass ( LTM ; Broadhurst et al. 2005; Zitrin et al. 2009, 2015) approach before the JWST ERO data release, and subsequently Pascale et al. (2022) present the LTM mass modeling with the JWST ERO data. Caminha et al. (2022) also develop the mass model of SMACS J0723 using Lenstool .  \nIn this paper, we construct an updated glafic (Oguri 2010, 2021) strong lens mass model of SMACS0723 using the new JWST ERO data. The magnification factors predicted by the updated glafic mass model are compared with those from the other existing mass models to evaluate the lens-model uncertainty. The glafic code performs the so-called parametric lens mass modeling,  \nwhere shapes of the mass distributions of the cluster are described by a superposition of a small number of lens mass components with known profile shapes, and parameters characterizing the lens mass components are determined so as to reproduce observed positions of multiple images.  \nAs a specific procedure, we largely follow the methodology described in Kawamata et al. (2016). We model the dark matter halo by an elliptical Navarro-FrenkWhite (NFW; Navarro et al. 1997) density profile with an approximation to speed up the calculation of its lensing property (see Oguri 2021). Model parameters associated with the NFW component are the mass, the center, the ellipticity and its position angle, and the concentration parameter. In addition to the main NFW halo, we place an additional NFW component whose center is fixed to a bright cluster member galaxy located at North-West, (R.A., Decl.)=(110 . 7928634, -73 . 4476417). We also fix the concentration parameter of the additional NFW component to c = 10, and fit its mass, ellipticity and position angle only. Cluster member galaxies selected by photometric redshifts from the RELICS HST data (Coe et al. 2019) are modeled by an elliptical pseudo Jaffe profile. In order to reduce the number of parameters, the velocity dispersion σ and the truncation radius r trunc of each cluster member galaxy are assumed to scale with its luminosity (in HST F814W band) as σ ∝ L 1 / 4 and r trunc ∝ L η with their normalization and η being treated as free parameters. In addition we include an external shear to improve the fitting. For multiple image sets without spectroscopic redshifts, we simultaneously fit their redshifts.  \nWe search for the best-fitting model by the standard χ 2 minimization, where χ 2 is computed from the differences between observed and model-predicted positions. We assume the positional error of 0 . '' 4, which is a typical positional accuracy achieved by the parametric strong lens mass modeling (Kawamata et al. 2016). The χ 2 is evaluated in the source plane, taking account of the full magnification tensor at each multiple image position (see Appendix 2 of Oguri 2010). Errors on model parameters are derived using the standard Markov Chain Monte Carlo technique. Multiple images are identified iteratively, starting with secure sets of multiple images that are obvious from their colors, morphologies, and redshifts, constructing a preliminary mass model with those sets of multiple images, and searching for new multiple sets with help of the preliminary mass model. In this work, we use conservative 12 sets of multiple images for our strong lens mass modeling, with the total number of multiple images of 38. These multiple image sets are mostly consistent with other work using different  \nTable 6. Summary of Magnification Factors Estimated by Various Mass Models  \nNote - (1) Name. (2) Spectroscopic or photometric redshift. (3)-(6) Magnification factors estimated by grafic , Mahler et al. (2022), Pascale et al. (2022), and Caminha et al. (2022).  \nlens modeling codes (Pascale et al. 2022; Mahler et al. 2022). We adopt spectroscopic redshifts for five sets of multiple images given in the literature (Golubchik et al. 2022; Pascale et al. 2022; Mahler et al. 2022). Our bestfitting model has χ 2 = 28 . 3 for degree of freedom of 32, representing a good fit. The root-mean-square (rms) of differences between observed and predicted multiple image positions is 0 . '' 35.  \nWith the updated glafic mass model, we calculate the magnification factors µ of our dropout galaxy candidates and the effective survey volume. Table 6 summarizes the magnification factors of our dropout galaxy candidate and spectroscopically confirmed galaxies at z > 7 calculated by glafic , Mahler et al. (2022), Pascale et al. (2022), and Caminha et al. (2022). We find that the magnification factor calculated by each model agree well typically within ∼ 20%.", '5.1. Sample Completeness': "To derive the rest-frame UV luminosity function, we estimate the completeness of our dropout galaxy selection in the same manner as previous studies (e.g., Ouchi et al. 2009; Ono et al. 2018; Harikane et al. 2022b). We conduct Monte Carlo simulations with real NIRCam images and artificial galaxies mocking high redshift galaxies. The mock high redshift galaxies follow the sizeM UV redshift distribution revealed with the HST legacy data sets for galaxies at z ∼ 0 -10 (Shibuya et al. 2015) that is extrapolated to our redshift ranges, where the sizeM UV distribution is the log-normal distribution. Our initial measurements of sizes for our z > 9 galaxy candidates are consistent with this assumption within the uncertainties. We adopt the S'ersic index n = 1 found in typical galaxies at z ∼ 5 -10 (Ono et al. 2013; Shibuya  \net al. 2015) and the flat distribution of the intrinsic ellipticity in the range of 0.0-0.8. Recent studies indicate that morphologies of z ∼ 9 -16 galaxies identified in the JWST datasets are consistent with these assumptions (e.g, Ono et al. 2022). The SEDs of the mock high redshift galaxies uniformly distribute over magnitude and redshift, and have a color distribution agreeing with the M UV -β UV relation observationally determined at z ∼ 8 (Bouwens et al. 2014), where β UV is the UV spectral slope index. The IGM absorption of Inoue et al. (2014) is applied to the SEDs, which produces absorption features in the wavelengths shorter than the Ly α line. We produce 100 artificial objects of the mock high redshift galaxies with IRAF mkobject in each redshift and magnitude bin, and place the artificial objects on the real JWST NIRCam images. With the images, we perform the object detections, photometry, the color selection, and the SED fitting in the same manner as Section 3. In the SMACS J0723 field, we consider the source magnification and multiply lensed images by using the mass model made with glafic described in Section 4. Finally we calculate the selection completeness as a function of magnitude and redshift, C ( m,z ), with the photometric catalogs of the artificial high redshift galaxies. Figure 11 presents examples of the selection completeness thus obtained. Although the average redshifts are z = 10 . 1 for F 115 W -dropouts, z = 13 . 8 for F 150 W -dropouts, and z = 18 . 7 for F 200 W -dropouts, we use the median of photometric redshifts of our selected candidates, z = 9 . 1 for F 115 W -dropouts, z = 12 . 0 for F 150 W -dropouts, and z = 16 . 3 for F 200 W -dropouts as the representative redshifts of our each dropout sample.  \nBased on the results of these selection completeness simulations, we estimate the survey volume per unit area as a function of apparent magnitude (Steidel et al. 1999),  \nV eff ( m ) = ∫ C ( m,z ) dV ( z ) dz dz, (10)  \nwhere dV ( z ) /dz is the differential comoving volume as a function of redshift. The space number density of our galaxy candidates that is corrected for incompleteness is calculated with the following equation:  \nψ ( m ) = n ( m ) V eff ( m ) , (11)  \nwhere n ( m ) is the surface number density of selected galaxies in an apparent magnitude bin of m . We convert the number density 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. Assuming a flat rest-frame UV continuum, we calculate the absolute UV  \nFigure 11. Selection completeness for our dropout galaxies. The purple, blue, and green curves show selection completeness for the F 115 W -, F 150 W -, and F 200 W -dropout galaxies whose rest-frame UV ( ∼ 1500 ˚ A) magnitudes are F 200 W = 27 . 0 mag, F 277 W = 27 . 0 mag, and F 356 W = 27 . 0 mag, respectively. Each selection window is smoothed by ∆ z = 1 . 0.  \n<!-- image -->  \nmagnitudes of galaxies from their apparent magnitudes in the bluest band not affected by the Lyman break, i.e., F 200 W , F 277 W , and F 356 W -bands for F 115 W , F 150 W , and F 200 W -dropout galaxy candidates, respectively. The 1 σ uncertainty is calculated by taking into account the Poisson confidence limit (Gehrels 1986) and the cosmic variance. We estimate the cosmic variance in the number densities using the bias values of z ∼ 7 galaxies obtained in Harikane et al. (2016), following the procedures in Somerville et al. (2004).", '5.2. Results': 'Figures 12 and 14 present our luminosity functions at z ∼ 9, 12, and 16 together with luminosity functions obtained by previous work including the latest JWST studies (Naidu et al. 2022b; Donnan et al. 2022; Finkelstein et al. 2022b; Bouwens et al. 2022b). Our measurements of the luminosity functions are summarized in Table 7. Comparing with the previous measurements of the luminosity functions, we find that our luminosity functions at z ∼ 9 and 12 agree well with those of the previous HST and JWST studies within the uncertainties, as shown in Figure 12. In Figure 14, we compare the luminosity function of our possible candidates at z ∼ 16 newly determined by this study with those available at lower redshifts at z ∼ 14 constrained by JWST. We confirm that these luminosity functions are comparable.  \nWe conduct χ 2 minimization fitting of the double power-law and Schechter functions to the luminosity functions that include the measurements at the bright  \n<!-- image -->  \n<!-- image -->  \nUV  \nFigure 12. UV luminosity functions at z ∼ 9 (left) and 12 (right). The red diamonds represent the number densities of our galaxy candidates, while the red arrows indicate the 1 σ upper limits. The errors include the cosmic variance (see text). The red solid and dashed lines are our best-fit double power-law and Schechter functions, respectively. In the left panel, the orange circles indicates the luminosity functions at z ∼ 9 obtained in Donnan et al. (2022) using JWST data, and the gray symbols and shades denote the results at z ∼ 9 derived by the previous studies using HST or ground-based telescope data, Oesch et al. (2013, squares), McLeod et al. (2016, right-pointing triangles), Morishita et al. (2018, pentagons), Stefanon et al. (2019, triangles), Bowler et al. (2020, circles), Bouwens et al. (2021, diamonds), Rojas-Ruiz et al. (2020, hexagons), Leethochawalit et al. (2022, down-pointing triangle), Finkelstein et al. (2022a, shade with dotted lines), and Bagley et al. (2022a, shade with dashed lines). In the right panel, the orange circles, the red circle, the orange down-pointing triangle, and the orange squares indicate the number density of galaxies at z ∼ 12, z ∼ 13, z ∼ 10 -13, and z ∼ 12 -13 reported by Donnan et al. (2022), Harikane et al. (2022a), Naidu et al. (2022b), and Bouwens et al. (2022b), respectively. The gray open symbols indicate the luminosity functions at z ∼ 10 obtained by McLeod et al. (2016, diamonds), Oesch et al. (2018, squares), Morishita et al. (2018, pentagons), and Bowler et al. (2020, circle). The green open star mark represents the number density of GN-z11 (Oesch et al. 2016). See Harikane et al. (2022a) for the estimate of the number density and the UV magnitude of GN-z11. Our estimated luminosity functions at z ∼ 9 and 12 agree well with previous HST and JWST results.Figure 13. Same as the right panel of Figure 12, but with the fitting results without the brightest datapoint in Harikane et al. (2022a).  \n<!-- image -->  \nend in the literature. In the fitting, we use the results of this study, Morishita et al. (2018), Bowler et al.  \n(2020), and Bouwens et al. (2021) for the z ∼ 9 luminosity function, the results of this study and Harikane et al. (2022a) for the z ∼ 12 luminosity function assuming that the UV luminosity function does not rapidly change at z ∼ 12 -13, and the result of this study for the z ∼ 16 luminosity function. We show the best-fit functions in Figures 12 and 14, and present the best-fit parameters in Table 8. At z ∼ 9, the χ 2 values of the fitting suggest that the double power-law function explains the luminosity functions ( χ 2 / dof = 2 . 3 / 9) better than the Schechter functions ( χ 2 / dof = 3 . 6 / 10), albeit with the moderately small difference of χ 2 ( ∼ 1 σ ). At z ∼ 12 and 16, we find no significant differences between the double power-law and Schechter functions in the χ 2 values, probably due to the large uncertainties of the measurements. At z ∼ 12, we also fit only the measurements of this study, excluding the brightest datapoint in Harikane et al. (2022a), as shown in Figure 13. The best-fit DPL and Schechter functions are slightly flatter than the fitting results with the datapoint in Harikane et al. (2022a) at the bright end.  \nFigure 14. UV luminosity function at z ∼ 16. The red diamond and the arrows represent the number density of our galaxy candidates and the 1 σ upper limits, respectively. For reference, we show the UV luminosity functions at the lower redshifts, z ∼ 14 (Donnan et al. 2022; orange filled circle), z ∼ 14 (Finkelstein et al. 2022b; orange filled square) z ∼ 12 (this study; gray open diamonds), z ∼ 13 (Harikane et al. 2022a; gray open circle), and z ∼ 10 -13 (Naidu et al. 2022b; gray open down-pointing triangle).  \n<!-- image -->  \nFigure 15 presents the redshift evolution of the luminosity function. We find the continuous decrease of luminosity functions from z ∼ 5 to z ∼ 12. We do not find a significant decrease from z ∼ 12 to 16 beyond the uncertainty. There is a hint of a small evolution from z ∼ 12 to 16, while the small number statistics do not allow us to conclude whether the evolutionary trend changes from z ∼ 5 -12 to 12 -16.  \nFigure 16 compares the observed luminosity functions at z ∼ 12 and 16 with those predicted by theoretical models (Dayal et al. 2014, 2019; Yung et al. 2020; Behroozi et al. 2020; Wilkins et al. 2022; Mason et al. 2022). At z ∼ 12, most of the models in Figure 16 explain the observational measurements in the faint magnitude range from -20 to -18 mag, while some models do not reproduce the moderately high number densities of the observational measurements at the bright magnitude of M UV < -20 mag. At z ∼ 16, most of the models cannot reproduce the observed number density of bright galaxies at M UV < -20 mag, except for the FLARES (Lovell et al. 2021; Vijayan et al. 2021; Wilkins et al. 2022) whose prediction at z ∼ 15 agrees with our number density estimate within uncertainties. Similarly, Figure 17 shows the predicted number of bright galaxies at z ∼ 12 -16 with M UV < -20 mag. Figure 17 indicates that the models underpredict the number of galaxies compared to the observation, although the sig-  \nFigure 15. Evolution of the UV luminosity functions from z ∼ 4 to z ∼ 16. The yellow, red, and pink diamonds represent our measurements of the luminosity functions at z ∼ 9, 12, and 16, respectively, whereas the red circle is the one obtained by Harikane et al. (2022a) at z ∼ 13. The orange, yellow, green, blue, purple, brown, and gray symbols indicate the luminosity functions at z ∼ 10, 9, 8, 7, 6, 5, and 4, respectively. The circles at z ∼ 4 -7 and 8 -10 are the data taken from Harikane et al. (2022b) and Bowler et al. (2020), respectively. The squares at z ∼ 4 -9 and z ∼ 10 are the data of Bouwens et al. (2021) and Oesch et al. (2018), respectively. The diamond at z ∼ 10 represents the result of McLeod et al. (2016). The lines denote the double powerlaw functions derived by the previous studies for z ∼ 4 -7 (Harikane et al. 2022b) and z ∼ 8 -13 (Bowler et al. 2020). For clarity, we shift the data point of Bowler et al. (2020) at z ∼ 10 by -0 . 2 mag.  \n<!-- image -->  \nUV  \nnificance is small and more data are needed to obtain the conclusion. This difference of the observations and models would suggest that the feedback effects in the models may be too strong to produce abundant bright galaxies, lower dust obscuration in these bright galaxies than the model assumptions, and/or that there exist hidden AGNs that produce radiation comparable with or more than stellar components of the galaxies (e.g., Bowler et al. 2014, 2020; Ono et al. 2018; Stevans et al. 2018; Shibuya et al. 2022; Harikane et al. 2022b; Pacucci et al. 2022; Mason et al. 2022), although there is also a possibility that this difference may be caused by other physical processes, as discussed in Section 6.', '5.3. Cosmic SFR Density': 'We derive the cosmic SFR densities at z ∼ 9, 12, and 16. We integrate the best-fit double power-law functions (Table 8) down to -17 mag, the same limit as previous studies (e.g., Bouwens et al. 2015; Oesch et al. 2018; Harikane et al. 2022a), and obtain the UV luminosity  \n<!-- image -->  \n<!-- image -->  \nUV  \nFigure 16. Comparison of the luminosity-function measurements with theoretical predictions and the empirical models at z ∼ 12 (left) and z ∼ 16 (right). The blue lines show the theoretical and empirical models obtained by Dayal et al. (2014, 2019, solid line), Yung et al. (2020, dotted line), Behroozi et al. (2020, dotted-dashed line), Wilkins et al. (2022, double-dotted dashed line), and Mason et al. (2022, dashed line; their no dust model). The red and orange symbols show observational results in the same manner as Figures 12 and 14. The red diamonds and arrows represent the measurements and upper limits obtained by this study. The orange circles, the red circle, the down-pointing orange triangle, and the orange square in the left (right) panel indicate the number densities reported by Donnan et al. (2022), Harikane et al. (2022a), Naidu et al. (2022b), and Bouwens et al. (2022b) (Finkelstein et al. (2022b)), respectively.Figure 17. Theoretical predictions for the number of bright galaxies at z ∼ 12 -16 with M UV < -20 mag detected in our survey area of ∼ 90 arcmin 2 . These numbers are based on the theoretical models of Dayal et al. (2014, 2019), Yung et al. (2020), Behroozi et al. (2020), Wilkins et al. (2022), and Mason et al. (2022). The red horizontal line with the shaded region indicates the number of observed galaxies at z ∼ 12 -16 with M UV < -20 mag ( N obs = 4 ± 2), which is higher than these model predictions.  \n<!-- image -->  \nTable 7. Obtained Luminosity Function at Each Redshift  \nNote -Errors and upper limits are 1 σ .  \nTable 8. Fit Parameters for Luminosity Functions  \nNote -Errors are 1 σ .  \nTable 9. Obtained Cosmic UV Luminosity Density and SFR Density  \nNote -Errors are 1 σ . ρ SFR , UV and ρ SFR are SFR densities without and with dust extinction correction, respectively.  \ndensities, ρ UV . We correct ρ UV for the dust extinction, following the attenuation-UV slope ( β UV ) relation (Meurer et al. 1999) and β UV -M UV relation at z = 8 in Bouwens et al. (2014). The choice of these assumptions (e.g., using the attenuation-UV slope law in de Barros et al. 2014 instead) does not affect our conclusions because the correction factor is very small ( glyph[lessorsimilar] 0 . 1 dex). We calculate SFRs from UV luminosities, L UV , corrected for dust extinction by the relation,  \nSFR ( M glyph[circledot] yr -1 ) = K UV L UV (erg s -1 Hz -1 ) , (12)  \nwhere K UV is the conversion factor that depends on the recent star-formation history, metal enrichment history, and the choice of the IMF. Here we apply K UV = 1 . 15 × 10 -28 M glyph[circledot] yr -1 / (erg s -1 Hz -1 ) that is used in Madau & Dickinson (2014). This value of K UV is valid for the Salpeter (1955) IMF, and consistent with the cosmic star-formation history and the evolved stellar metallicity (10 -0 . 15 z Z glyph[circledot] ; Madau & Dickinson 2014) up to z ∼ 10. Table 9 summarizes our measurements of the cosmic UV  \nluminosity density, SFR densities without and with dust extinction correction at each redshift.  \nFigure 18 presents the cosmic SFR density evolution. In this figure, we show the cosmic SFR density measurements at z ∼ 0 -10 obtained by previous studies, all of which are converted to the calibration of Madau & Dickinson (2014) with the Salpeter (1955) IMF (Equation (12)). We confirm that our SFR density at z ∼ 9 is consistent with the previous measurements. We compare the observational measurements of the SFR densities with the constant star-formation efficiency ( SFR/ ˙ M h ( z ) = const . ) model (Harikane et al. 2022b) together with the extrapolation of the Madau & Dickinson (2014) estimates at z = 0 -8. We find that the cosmic SFR densities significantly decrease from z ∼ 9 to 12. A decrease of the cosmic SFR densities may exist from z ∼ 12 to 16, while the decrease is not larger than the errors. Interestingly, The constant star-formation efficiency model explains the evolution of the cosmic SFR densities up to z ∼ 10 (Harikane et al. 2022b), while  \nFigure 18. Cosmic SFR density evolution. The red circles represent the cosmic SFR densities obtained by our study, with the double power-law luminosity functions integrated down to M UV = -17 mag. The black circles indicate the cosmic SFR densities derived by Madau & Dickinson (2014), Finkelstein et al. (2015b), McLeod et al. (2016), Bhatawdekar et al. (2019), and Bouwens et al. (2020). The orange circles are results in Donnan et al. (2022). The blue dashed curve is the best-fit function of the cosmic SFR densities in Harikane et al. (2022b, their Equation 60). In Harikane et al. (2022b), they assume the constant star formation efficiency at z > 10, resulting in the power-law decline with ρ SFR ∝ 10 -0 . 5(1+ z ) . The gray dashed curve shows the best-fit function at z glyph[lessorsimilar] 8 determined by Madau & Dickinson (2014) extrapolated to z > 8. All results are converted to those of the Salpeter (1955) IMF.  \n<!-- image -->  \nour measurement at z ∼ 12 is higher than the model prediction beyond the uncertainty. Moreover, there is a hint of a high cosmic SFR density at z ∼ 16 above the model, although it is not statistically significant due to the large error. Such higher SFR densities than the constant efficiency model at z ∼ 15 is actually consistent with observations of Balmer break galaxy candidates at z ∼ 6 (Mawatari et al. 2020b).', '6.1. Possible High Cosmic SFR Density at z > 10': 'Our observational measurements suggest that the SFR densities at z ∼ 12 -16 are higher than the constant starformation efficiency model of Harikane et al. (2022b). Although the constant star-formation efficiency model well explains the cosmic SFR densities at z ∼ 0 -10, this model underpredicts those at z ∼ 12 -16. Here we discuss the following three possibilities that explain the observed high SFR densities at z ∼ 12 -16.  \n- (A) No star formation suppression at the prereionization epoch. The universe at z ∼ 12 -16 is at the pre-reionization epoch when the IGM is highly neutral (Ouchi et al. 2020; Robertson 2021). At the epoch of reionizaton (EoR; z ∼ 6 -12) and the epoch of post-reionization (postEoR; z glyph[lessorsimilar] 6), galaxies and AGN produce UV radiation by their star-formation and nuclear activity, and produce strong UV background radiation. The UV background radiation heats up Hi gas in low-mass halos of M h glyph[lessorsimilar] 10 8 -9 M glyph[circledot] with negligible Hi self-shielding, suppressing star-formation at the EoR and post-EoR (Barkana & Loeb 2000; Susa & Umemura 2004; Hoeft et al. 2006; Pawlik & Schaye 2009; Mesinger et al. 2009; Sawala et al. 2010; Bland-Hawthorn et al. 2015). Although the halo masses of our galaxies at z ∼ 12 -16 are unknown, the maximum halo mass existing at z ∼ 15 given the survey volume of this study is M h glyph[similarequal] 3 × 10 9 M glyph[circledot] in the structure formation model with the Planck cosmology (Behroozi et al. 2020). In other words, most of halos (with glyph[lessorsimilar] 10 9 M glyph[circledot] ) at z ∼ 15 are not affected by the UV background at the pre-reionization epoch, while the similar halos at z glyph[lessorsimilar] 10 experience the suppression of star formation by the UV background at the EoR and postEoR. To test whether this effect can quantitatively explains the observed SFR densities, we construct a model of the SFR density evolution including the enhancement of the star formation that is free from the suppression by the UV background at the pre-EoR. We use a model in Barkana & Loeb (2000) with a reionization redshift of z reion = 13, and multiply the prediction of the constant star formation efficiency model (Harikane et al. 2022b) by a factor of the star formation rate enhancement due to no suppression by the UV background in Barkana & Loeb (2000). The left panel of Figure 19 presents this hybrid model including the effect of star formation enhancement at pre-EoR, which reproduces the observed SFR densities at z ∼ 12 -16 within uncertainties. This agreement indicates a possibility that the star formation efficiency at z ∼ 12 -16 is higher than those at z glyph[lessorsimilar] 10 due to no suppression of the star formation activity at the pre-EoR.\n- (B) Presence of AGN activity. Another possibility is that a large fraction of the observed UV luminosity densities at z ∼ 12 -16 is produced by AGN, and there is no excessive SFR densities at z ∼ 12 -16 beyond the constant star-formation efficiency model. This is an interesting scenario that  \n<!-- image -->  \nFigure 19. Possible scenarios to explain the observed SFR densities at z > 10. (Left:) Scenario of no star formation suppression at pre-reionization epoch. At the reionization epoch and after that, star formation in low-mass halos is suppressed by strong UV background radiation, while before the reionization epoch such a suppression of star formation activity does not occur. The upper edge of the blue shaded region indicates the enhancement of the star formation by this effect (Barkana & Loeb 2000), which explains the observed SFR densities (see texts for details). (Right:) Scenario of Pop III star formation. As shown in Figure 20, Pop III stellar populations with a top-heavy IMF produces a significant amount of UV photons at a given SFR, resulting in the overestimates of the SFR densities if we use the canonical UV-SFR conversion factor. The red and orange filled circles at z > 10 . 5 are SFR densities calculated based on the conversion factor for a Pop III stellar population with a top-heavy IMF (the PopIII.1 model in Figure 20), which agree well with the constant star formation efficiency model in Harikane et al. (2022b). The open circles are SFR densities based on the canonical conversion factor.  \n<!-- image -->  \nmitigates the existence of supermassive black holes (SMBH) at z ∼ 7 (Mortlock et al. 2011, Ba˜nados et al. 2018, Wang et al. 2021) by efficient gas accretion on SMBHs creating AGNs, while a standard gas accretion limited by the Eddington accretion rate does not explain the existence of the SMBHs at z ∼ 7. However, our z ∼ 12 -16 candidates except for GL-z12-1 show extended morphologies (Section 3.5). Thus the fraction of AGN radiation dominated galaxies is as small as ∼ 10% (= 1 / 10) at z ∼ 12 -16. Although the excessive SFR density estimate at z ∼ 16 is unclear due to the small statistics, the one at z ∼ 12 cannot be explained by AGN activity.  \n(C) A top-heavy IMF. The third possibility is an overestimate of the SFR density due to a top-heavy IMF possibly with the Population III (Pop III) stellar population. In our estimate of the SFR density, we use the canonical UV luminosityto-SFR conversion factor of K UV = 1 . 15 × 10 -28 M glyph[circledot] yr -1 / (erg s -1 Hz -1 ), which is for the Salpeter (1955) IMF, while K UV depends on starformation history, metallicities, and IMFs (e.g., Madau & Dickinson 2014; Tacchella et al. 2018). Indeed in the early universe, the IMF is expected to be more top-heavy because of a lower metallic-  \nty or a higher CMB temperature (e.g., Omukai et al. 2005; Chon et al. 2022), resulting in a higher Jeans mass, especially for Pop III stellar populations (e.g., Hirano et al. 2014, 2015). To test whether this effect can explain the observed densities, we calculate the UV-to-SFR conversion factor, K UV , for different metallicity and IMF assumptions using Yggdrasil (Zackrisson et al. 2011). Figure 20 presents K UV for different metallicities and IMFs as a function of stellar age. We find that Pop III stellar populations with top-heavy IMFs (PopIII.1 and PopIII.2 in Yggdrasil ) produce ∼ 3 -4 times more UV photons than the canonical assumption given the SFR, because nebular continuum emission boosts the UV luminosity as discussed in previous studies (e.g., Zackrisson et al. 2008; Schaerer & de Barros 2009, 2010). This low conversion factor reduces the SFR density estimates at z ∼ 12 -16 as shown in the right panel of Figure 19, resulting in the SFR densities consistent with the constant star formation efficiency model.  \nBased on these discussions, we conclude that (A) no star formation suppression at pre-reionization epoch or (C) a top-heavy IMF with a Pop III-like star formation can explain the observed high SFR densities at  \nFigure 20. UV luminosity-SFR conversion factor, K UV , for various metallicities as a function of stellar age. The green, cyan, and blue curves show the conversion factor for metallicities of Z = 0 . 02, 0 . 004, 0 . 0004, respectively, calculated with Yggdrasil (Zackrisson et al. 2011) assuming a constant star formation history with a unity gas covering fraction. These factors are values for a UV luminosity at 1500 ˚ A in the Salpeter (1955) IMF in the interval of 0 . 1 -100 M glyph[circledot] . Note that the original outputs from Yggdrasil are for the Kroupa (2001) IMF, and we correct for the IMF difference by multiplying the outputs by 1.49. The solid (PopIII.1) and dashed (PopIII.2) purple curves show the conversion factors for Pop III stellar populations with an extremely top-heavy IMF (50-500 M glyph[circledot] , the Salpeter 1955 slope) and a moderately top-heavy IMF (log-normal with characteristic mass of M c = 10 M glyph[circledot] , dispersion σ = 1 M glyph[circledot] , and wings extending from 1 -500 M glyph[circledot] ). If galaxies at z > 10 are dominated by Pop III stellar populations, the conversion factor is significantly lower than the typically assumed value ( K UV = 1 . 15 × 10 -28 M glyph[circledot] yr -1 / (erg s -1 Hz -1 ); the black line), resulting in the overestimate of the SFR.  \n<!-- image -->  \nz ∼ 12 -16. These possibilities can be further investigated by follow-up observations with JWST/MIRI covering a longer wavelength than the Balmer break to obtain the robust stellar mass measurements and star formation history, or with JWST/NIRSpec and MIRI spectroscopy to search for signatures of Pop III-like stellar populations and AGN activity.', '6.2. Properties of Luminous Galaxy Candidates': 'In this study, we have found several luminous galaxy candidates at the early epoch of z ∼ 10 -16, when the age of the universe is only ∼ 200 -500 Myrs after the Big Bang. Here we discuss physical properties of these luminous galaxy candidates.  \nTable 10. SFRs and Stellar Masses of Luminous Galaxy Candidates with M UV < -19 . 5 mag  \nNote -Assuming the Chabrier (2003) IMF and metallicity of Z = 0 . 2 Z glyph[circledot] . The SFR is averaged over the past 50 Myr in the same manner as Tacchella et al. (2022). See Section 3.3 for the details of the SED fitting.  \n- † This candidate shows a compact morphology indicative of AGN activity, while the profile is spatially extended more than the PSF (Section 3.5).  \nTable 10 summarizes SFRs and stellar masses of six galaxy candidates whose UV magnitudes are brighter than M UV = -19 . 5 mag, constrained by the SED fitting in Section 3.3 assuming the Chabrier (2003) IMF. Our estimates of SFRs and stellar masses agree with previous estimates by Naidu et al. (2022b), Donnan et al. (2022), and Finkelstein et al. (2022b), indicating that these luminous galaxies are very massive with stellar masses as high as M ∗ ∼ (1 -10) × 10 8 M glyph[circledot] at z ∼ 10 -16. While the contributions from AGN radiation to the SEDs may be suspected in one of the objects, GL-z12-1 (see Sections 3.5 and 6.1), at least the rest of the objects (i.e. ∼ 80% of the bright z ∼ 10 -16 galaxies) would be truly stellar massive. Although the NIRCam photometry is limited to < 5 µ m and does not trace the SEDs beyond the Balmer break (4000 ˚ A) corresponding to 5 -7 µ m in the observed frame at z ∼ 10 -16, these stellar mass estimates provide rough lower limits that miss the contribution from old stellar populations beyond the Balmer break, given high specific SFRs of these galaxy candidates, SFR/M ∗ ∼ 10 -8 yr -1 .  \nHere is a question how these galaxies with the large stellar masses form at this early epoch of z ∼ 10 -16. To discuss the formation scenario of these massive galaxy candidates, we estimate the stellar-to-halo mass ratio (SHMR) of these galaxies. Using the abundance matching technique in the same manner as Harikane et al. (2016, their Equation (66)), we estimate the halo mass of the most massive halo that can be observed with the survey volume in this study, resulting in 5 × 10 10 M glyph[circledot] and 5 × 10 9 M glyph[circledot] at z ∼ 12 and z ∼ 16, respectively.  \nFigure 21. Stellar masses of our galaxy candidates as a function of redshift. The red filled diamonds show the stellar mass estimates for six luminous galaxy candidates with M UV < -19 . 5 mag at z ∼ 10 -16 (Table 10), and open red circles are results for other candidates. The gray shaded region indicates the stellar mass whose number density is below the observational limit, calculated from the cosmic baryon fraction (Ω b / Ω m = 0 . 16) and the maximum halo mass that can be observed with the survey volume of this study. The black and blue curves indicate the stellar masses calculated from the maximum halo mass with M ∗ /M h = 0 . 1 and the maximum M ∗ /M h value at each redshift in Behroozi et al. (2020), respectively. The massive stellar masses ( M ∗ ∼ 10 9 M glyph[circledot] ) of the two z ∼ 16 candidates can be explained by a very high SHMR of M ∗ /M h = 0 . 1, indicating a star formation efficiency as high as ∼ 60%. The other four luminous candidates at z ∼ 10 -13 also show higher stellar masses compared to the predictions from Behroozi et al. (2020).  \n<!-- image -->  \nFrom the stellar mass estimates discussed above, the SHMRs of z ∼ 12 and z ∼ 16 galaxies are ∼ 0 . 01 and ∼ 0 . 1, respectively. Because the cosmic baryon to dark matter density ratio is Ω b / Ω m = 0 . 16 (Planck Collaboration et al. 2020), the SHMRs of the z ∼ 16 galaxies reach ∼ 60% of the cosmic baryon fraction, as shown in Figure 21. In other words, more than a half of baryon gas in the halos is converted to stars, which is unlikely found in lower-redshift and present-day galaxies whose SHMRs are ∼ 0 . 02 -0 . 03 at maximum (e.g., Harikane et al. 2016; Behroozi et al. 2019). A similar conclusion is obtained from the comparison of the UV luminosity functions (Inayoshi et al. 2022). However, theoretical models predict such efficient star formation at the pre-reionization epoch where 70 -80% of baryon are converted stars (Susa & Umemura 2004) in halos with 10 8 -9 M glyph[circledot] masses in a few hundred Myr, when the UV background radiation is too weak to suppress star forma-  \non (A) in Section 6.1). The other four galaxies at z ∼ 11 -14 also show higher stellar masses compared to the predictions from the maximum SHMR in Behroozi et al. (2020), indicating elevated star formation efficiencies, probably due to no suppression of star formation activity at the pre-reionization epoch. Another possibility is that the SFR and the stellar mass of these bright galaxies are overestimated due to the assumption of the IMF and metallicity in the SED fitting, as discussed in Section 6.1 (discussion (C)). Indeed, if we assume that the stellar population of these galaxies are dominated by Pop III with a top-heavy IMF, the SFR and stellar mass are reduced by a factor of ∼ 3 -4, more comparable to the observed SHMRs at lower redshifts. These comparisons, together with the discussions in Section 6.1, indicate that the observed properties of z ∼ 10 -16 galaxies (i.e., high cosmic SFR densities and massive stellar masses) can be explained by either no star formation suppression by UV background radiation at the pre-reionization epoch or a top heavy IMF possibly with a Pop III-like stellar population.', '7. SUMMARY': "In this paper, we have conducted comprehensive analyses for the JWST/NIRCam images taken by the JWST ERO SMACS J0723, Stephan's Quintet, ERS GLASS, and CEERS projects, covering a total of ∼ 88 . 7 arcmin 2 , in conjunction with the supports of the ERO SMACS J0723 NIRSpec spectra. We reduced the NIRCam datasets using the new calibration parameters released in October, based on calibration observations of three different standard stars placed in all of the 10 NIRCam detectors. Our major findings are summarized below:  \n- 1. We have selected dropout galaxy candidates at z ∼ 9, z ∼ 12, and z ∼ 16 showing significant continuum breaks in the NIRCam F 115 W , F 150 W , and F 200 W -bands, respectively, by the color criteria, confirming clear non-detections in the band(s) whose wavelength is shorter than the continuum breaks including the F 090 W band (Section 3.2, Figure 2). Because we have found that a weak photoz criterion of ∆ χ 2 > 4 cannot remove a number of foreground interlopers on the bases of the JWST simulation data produced by the CEERS project team (Figure 3), we apply a stringent photoz determination criterion of ∆ χ 2 > 9 with the prospector code for our galaxy selection. We thus identify 13, 8, and 2 dropout galaxy candidates at z ∼ 9, z ∼ 12, and z ∼ 16, respectively (Table 2). We confirm that our photometric redshifts agree well with the spectroscopic redshifts, by applying our photometric redshift tech-  \nat z spec ∼ 8 -9 found by the ERO NIRSpec observations (Figures 4 and 5).  \n- 2. We have thoroughly compared our dropout galaxy candidates with other high redshift galaxies reported in a number of recent studies in the ERO SMACS J0723 and the ERS GLASS+CEERS NIRCam fields. We have summarized the candidates so far claimed in the literature together with our dropouts in Tables 3-5. For bright galaxy candidates, we find that a reasonable fraction of galaxies are commonly selected in our and previous studies. We confirm that, among all of the candidates, our dropout galaxies show the significant Ly α continuum breaks and flat UV continua with non-detections of continua whose wavelengths are shorter than the break (Figures 7-8), and conclude that we do not miss many reliable candidates in the redshift range of z ∼ 9 -16 in our selection.\n- 3. We have derived the UV luminosity functions at z ∼ 9, 12, and 16 (Figures 12 and 14). The UV luminosity functions at z ∼ 9 and 12 agree with those of previous HST and JWST studies within uncertainties including the cosmic variance, and the UV luminosity function at z ∼ 16 is newly constrained. The double power-law function is preferred to the Schechter function at z ∼ 9, albeit with the moderately small difference of χ 2 .\n- 4. The cosmic SFR densities at z ∼ 9, 12, and 16 are derived by the integration of the best-fit UV luminosity functions (Figure 18). By the comparisons with the previous low-redshift determinations of cosmic SFR densities, we find that the cosmic SFR densities significantly decrease from z ∼ 9 to 12. A decrease of the cosmic SFR densities may exist from z ∼ 12 to 16, while the decrease is not larger than the errors. Our measurements of the cosmic SFR density at z ∼ 12 is higher than predictions from the constant star-formation efficiency model (Harikane et al. 2022b), while the model explains the cosmic star-formation history at z glyph[lessorsimilar] 10. Moreover, there is a hint of a high cosmic SFR density at z ∼ 16 above the model, although it is not statistically significant due to the large error.\n- 5. There are several luminous and massive galaxy candidates with M UV < -19 . 5 mag at the early epoch of z ∼ 10 -16, when the age of the universe is only ∼ 200 -500 Myrs after the Big Bang (Figure 21). We confirm that our stellar mass estimates are comparable with those of the previous studies. Although one of the objects  \nmay have contributions of UV radiation from an AGN suggested by their morphologies, a majority ( ∼ 80%) of the galaxies may be truly stellar massive. By the comparisons with the structure formation models that provide the upper limits of the dark-matter halo masses observed in this study, the SHMR of the luminous galaxy candidates at z ∼ 16 is M ∗ /M h ∼ 0 . 1, corresponding to ∼ 60% of the baryon to dark matter density ratio in the Planck cosmology, indicating that most of baryon may be converted to stars, unlike lower-redshift and present-day galaxies with a reasonably small SHMR up to M ∗ /M h glyph[similarequal] 0 . 02 -0 . 03 (e.g., Harikane et al. 2016; Behroozi et al. 2019). The other candidates at z ∼ 10 -13 also have stellar masses more massive than predictions from the maximum SHMR in Behroozi et al. (2020).  \n- 6. This study identifies two interesting observational properties of galaxies at z ∼ 10 -16, the cosmic SFR densities higher than the constant starformation efficiency model and the existence of the UV-luminous galaxies with high stellar masses. The possibility of the AGN contribution can be ruled out, because the small fraction of galaxies have compact morphologies suggesting no dominant radiation from the AGN activity. Instead, there are two scenarios that explain the observational properties (Figure 19). One scenario is that the UV background radiation does not suppress the star formation at the pre-reionization epoch unlike at the EoR and post-EoR. Efficient star formation may take place at z ∼ 10 -16, producing the high cosmic SFR densities and the stellar massive galaxies. The other scenario is that a top-heavy IMF possibly with Pop III (or similarly metal poor) stellar populations produces strong UV radiation. The strong UV radiation may result in the overestimates of SFR densities above the constant star-formation efficiency model and of the stellar mass of the luminous galaxies. Further observational and theoretical studies are needed to test these two scenarios.", 'ACKNOWLEDGMENTS': 'We thank the anonymous referee for a careful reading and valuable comments that improved the clarity of the paper. We thank the CEERS team, especially Micaela Bagley and Steven Finkelstein, for providing many useful scripts for the NIRCam data reduction and datasets of the CEERS simulated images, and L. Y. Aaron Yung for providing SEDs of the mock galaxies. We are grateful to Rychard Bouwens, Seiji Fujimoto, Kohei Inayoshi, Akio Inoue, Tadayuki Kodama, Sandro Tacchella, Kenichi Tadaki, and Masayuki Umemura for useful comments and discussions. We thank Pratika Dayal and L. Y. Aaron Yung for sharing their data of the luminosity function. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with programs 2732, 2736, 1324, and 1345. The authors acknowledge the ERO, GLASS, and CEERS teams led by Klaus M. Pontoppidan, Tommaso Treu, and Steven L. Finkelstein, respectively, for developing their observing programs with a zero-exclusive-access period. This publication is based upon work supported by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and KAKENHI (20H00180, 20H00181, 20H05856, 21K13953, 21H04467, 22H01260) through Japan Society for the Promotion of Science. This work was supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo.  \nSoftware: PANHIT (Mawatari et al. 2020a), Prospector (Johnson et al. 2021), SExtractor (Bertin & Arnouts 1996), SWarp (Bertin et al. 2002), Yggdrasil (Zackrisson et al. 2011)'}