We present new measurements of restUV luminosity functions and angular correlation functions from 4100221 galaxies at z  27 identified in the SubaruHyper SuprimeCam survey and CFHT Large Area Uband Survey. The obtained luminosity functions at z  47 cover a very wide UV luminosity range of  0.002mbox2000LmathrmUV  combined with previous studies confirming that the dropout luminosity function is a superposition of the active galactic nucleus AGN luminosity function dominant at M SUBUVSUB  24 mag and the galaxy luminosity function dominant at M SUBUVSUB  22 mag consistent with galaxy fractions based on 1037 spectroscopically identified sources. Galaxy luminosity functions estimated from the spectroscopic galaxy fractions show the brightend excess beyond the Schechter function at 2 levels possibly made by inefficient mass quenching low dust obscuration andor hidden AGN activity. By analyzing the correlation functions at z  26 with HOD models we find a weak redshift evolution within 0.3 dex of the ratio of the star formation rate SFR to the dark matter accretion rate mathrmSFRdotMrmh  indicating the almost constant star formation efficiency at z  26 as suggested by our earlier work at z  47. Meanwhile the ratio gradually increases with decreasing redshift at z lt 5 within 0.3 dex which quantitatively reproduces the cosmic SFR density evolution suggesting that the redshift evolution is primarily driven by the increase of the halo number density due to the structure formation and the decrease of the accretion rate due to the cosmic expansion. Extrapolating this calculation to higher redshifts assuming the constant efficiency suggests a rapid decrease of the SFR density at z gt 10 with  10SUP0.51zSUP which will be directly tested with the James Webb Space Telescope.
['Galaxy evolution', 'Galaxy formation', 'High-redshift galaxies', '594', '595', '734', 'Astrophysics - Astrophysics of Galaxies']
GOLDRUSH. IV. Luminosity Functions and Clustering Revealed with  4000000 Galaxies at z   27 GalaxyAGN Transition Star Formation Efficiency and Implication for Evolution at z gt 10
{'GOLDRUSH. IV. Luminosity Functions and Clustering Revealed with ∼ 4,000,000 Galaxies at z ∼ 2 -7 : Galaxy-AGN Transition, Star Formation Efficiency, and Implication for Evolution at z > 10': "Yuichi Harikane 1,2 , Yoshiaki Ono 1 , Masami Ouchi 3,4,1,5 , Chengze Liu 6 , Marcin Sawicki 7,8 , Takatoshi Shibuya 9 , Peter S. Behroozi 10 , Wanqiu He 3 , Kazuhiro Shimasaku 11,12 , Stephane Arnouts 13 , Jean Coupon 14 , Seiji Fujimoto 15,16 , Stephen Gwyn 17 , Jiasheng Huang 18,19 , Akio K. Inoue 20,21 , Nobunari Kashikawa 11,12 , Yutaka Komiyama 3,4 , Yoshiki Matsuoka 22 , and Chris J. Willott 17  \n- 1 Institute for Cosmic Ray Research, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8582, Japan 2 Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK 3  \nNational Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  \n- 4 Graduate University for Advanced Studies (SOKENDAI), 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan\n- 5 Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), The University of Tokyo, 5-1-5 Kashiwanoha,Kashiwa, Chiba, 277-8583, Japan\n- 6 Department of Astronomy, School of Physics and Astronomy, and Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai 200240, P. R. China\n- 7 Institute for Computational Astrophysics and Department of Astronomy and Physics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada  \n8 Canada Research Chair  \n9 Kitami Institute of Technology, 165, Koen-cho, Kitami, Hokkaido 090-8507, Japan  \n- 10 Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  \nDepartment of Astronomy, School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  \n12  \nResearch Center for the Early Universe, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  \n13  \n14  \nAix Marseille University, CNRS, CNES, Laboratoire d'Astrophysique de Marseille, Marseille, France  \nAstronomy Department, University of Geneva, Chemin d'Ecogia 16, CH-1290 Versoix, Switzerland  \n15 Cosmic Dawn Center (DAWN), Copenhagen Denmark  \n16 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, Copenhagen, Denmark  \nNRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC, V9E 2E7, Canada  \n18 CASSACA, National Astronomical Observatories of China  \n19 Harvard-Smithsonian Centre for Astrophysics, USA  \n- 20 Department of Physics, School of Advanced Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555\n- 21 Waseda Research Institute for Science and Engineering, Faculty of Science and Engineering, Waseda University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan and\n- 22 Research Center for Space and Cosmic Evolution, Ehime University, Bunkyo-cho, Matsuyama, Ehime 790-8577, Japan Accepted for Publication in ApJS", 'Abstract': 'We present new measurements of rest-UV luminosity functions and angular correlation functions from 4,100,221 galaxies at z ∼ 2 -7 identified in the Subaru/Hyper Suprime-Cam survey and CFHT Large-Area U -band Survey. The obtained luminosity functions at z ∼ 4 -7 cover a very wide UV luminosity range of ∼ 0 . 002 -2000 L ∗ UV combined with previous studies, confirming that the dropout luminosity function is a superposition of the AGN luminosity function dominant at M UV glyph[lessorsimilar] -24 mag and the galaxy luminosity function dominant at M UV glyph[greaterorsimilar] -22 mag, consistent with galaxy fractions based on 1037 spectroscopically-identified sources. Galaxy luminosity functions estimated from the spectroscopic galaxy fractions show the bright end excess beyond the Schechter function at glyph[greaterorsimilar] 2 σ levels, possibly made by inefficient mass quenching, low dust obscuration, and/or hidden AGN activity. By analyzing the correlation functions at z ∼ 2 -6 with halo occupation distribution models, we find a weak redshift evolution (within 0 . 3 dex) of the ratio of the star formation rate (SFR) to the dark matter accretion rate, SFR/ ˙ M h , indicating the almost constant star formation efficiency at z ∼ 2 -6, as suggested by our earlier work at z ∼ 4 -7. Meanwhile, the ratio gradually increases with decreasing redshift at z < 5 within 0 . 3 dex, which quantitatively reproduces the cosmic SFR density evolution, suggesting that the redshift evolution is primarily driven by the increase of the halo number density due to the structure formation, and the decrease of the accretion rate due to the cosmic expansion. Extrapolating this calculation to higher redshifts assuming the constant efficiency suggests a rapid decrease of the SFR density at z > 10 with ∝ 10 -0 . 5(1+ z ) , which will be directly tested with JWST .  \nKey words: galaxies: formation - galaxies: evolution - galaxies: high-redshift', '1. Introduction': "Studying statistical properties of galaxies is important to understand the overall picture of galaxy formation and evolution. To quantify galaxy build-up in the early universe, many studies have investigated luminosity functions (i.e., one-point statistics) and angular correla-  \ntion functions (i.e., two-point statistics) of high redshift galaxies. The luminosity function represents the volume density of galaxies as a function of the luminosity. Since galaxies form in dark matter halos, the luminosity function is related to the dark matter mass function and baryonic physics of galaxy formation. Studying the shape and evolution of the luminosity function in the high redshift universe allows us to obtain key insights into the star  \n11  \n17  \nformation and feedback processes.  \nGreat progress has been made in determining luminosity functions of high redshift galaxies, especially in the rest-frame ultraviolet (UV), which is redshifted to the optical wavelength at z ∼ 4 -7 easily accessible from ground-based telescopes. Since the time-averaged unobscured star formation rate (SFR) of galaxies is proportional to the luminosity of galaxies in the rest-frame UV, the UV luminosity function provides us with a measure of how quickly galaxies grow with cosmic time. Analyses of galaxies in deep blank fields including the Hubble Ultra Deep field (HUDF) have resulted in identifying ∼ 20 , 000 galaxy candidates at z ∼ 2 -10 down to the absolute UV magnitude of M UV ∼ -17 mag (e.g., Oesch et al. 2010, Ellis et al. 2013; Schenker et al. 2013; McLure et al. 2013; Bouwens et al. 2015, 2019, 2021; Finkelstein et al. 2015b; Parsa et al. 2016; Mehta et al. 2017). In addition, the gravitational lensing by galaxy clusters has allowed us to probe even fainter galaxies and constrain the faint-end slope of the UV luminosity function (e.g., Ishigaki et al. 2015, 2018; Oesch et al. 2015, 2018, Atek et al. 2015, 2018; Castellano et al. 2016; Alavi et al. 2016; McLeod et al. 2016, Bouwens et al. 2017), although the impact of magnification uncertainties should be correctly considered (see Bouwens et al. 2017, Atek et al. 2018).  \nInvestigating the bright-end of the luminosity function is also important. Previously the luminosity function is thought to follow the Schechter function (Schechter 1976), which is derived from the shape of the halo mass function (Press & Schechter 1974) with several modifications. The Schechter function has an exponential cutoff at the bright end, which is possibly attributed to several different mechanisms such as heating from an active galactic nucleus (AGN; e.g., Binney 2004; Scannapieco & Oh 2004; Granato et al. 2004; Croton et al. 2006; Bower et al. 2006), inefficiency of gas cooling in massive dark matter haloes due to virial shock heating (e.g., Binney 1977; Rees & Ostriker 1977; Silk 1977; Benson et al. 2003), and dust obscuration which becomes substantial for the most luminous galaxies (e.g., Wang & Heckman 1996; Adelberger & Steidel 2000; Martin et al. 2005, Bowler et al. 2020). However, recent studies based on wide area surveys have reported an overabundance of objects at the bright end of UV luminosity functions beyond the Schechter function (bright end excess, e.g., Bowler et al. 2014, 2015, 2017, 2020; Stefanon et al. 2017, 2019; Ono et al. 2018; Stevans et al. 2018; Morishita et al. 2018; Adams et al. 2020; Moutard et al. 2020; Finkelstein et al. 2021). These studies suggest that the bright-end ( ≤ -23 mag) of the luminosity function is contributed by faint quasars or AGNs, at least at z ∼ 4 -7 (Ono et al. 2018; Stevans et al. 2018; Adams et al. 2020). In addition, Ono et al. (2018) calculate the galaxy UV luminosity function that is estimated by the subtraction of the AGN contribution, and report that the galaxy luminosity function still shows a bright end excess beyond the Schechter function at z ∼ 4 -7. They claim that this bright end excess implies inefficient mass quenching (e.g., the AGN feedback, virial shock heating) in these high redshift galaxies, or significant number of merging or gravitationally-lensed galaxies at the bright end (see also; e.g., Bowler et al. 2014; Bouwens et al. 2015).  \nTogether with studying luminosity functions, the clustering analysis with the angular correlation function is  \nimportant to understand the connection between galaxies and their dark matter halos. The galaxy-dark matter halo connection is investigated with the weak lensing analysis (e.g., Leauthaud et al. 2012; More et al. 2015; Coupon et al. 2015), the abundance matching/empirical model (e.g., Behroozi et al. 2013, 2019; Moster et al. 2013, 2018; Finkelstein et al. 2015a), and the clustering analysis (e.g., Harikane et al. 2016, 2018a; Ishikawa et al. 2017, 2020; Cowley et al. 2018; Qiu et al. 2018; Cheema et al. 2020). Since the weak lensing analysis cannot be applied at z glyph[greaterorsimilar] 2 due to the limited number of the background galaxies and their lower image quality with current observational datasets (but see also Miyatake et al. 2021), the clustering analysis is a crucial tool for estimating the dark matter halo mass of high redshift galaxies. Many studies have investigated the dark matter halos of high redshift galaxies as a function of their redshifts and UV luminosities (e.g., Ouchi et al. 2004, 2005; Hildebrandt et al. 2009; Savoy et al. 2011; Bian et al. 2013; BaroneNugent et al. 2014; Harikane et al. 2016, 2018b; Hatfield et al. 2018). These studies reveal that the more UV luminous galaxies reside in more massive halos.  \nRecently, Harikane et al. (2018b) have identified a tight relation between the ratio of the SFR to the dark matter accretion rate, SFR/ ˙ M h , and the halo mass, M h , over z ∼ 4 -7, suggesting the existence of a fundamental relation between the galaxy growth and its dark matter halo assembly. This redshift-independent relation indicates that the star formation efficiency does not significantly change at z ∼ 4 -7, and star formation activities are regulated by the dark matter mass assembly. Several studies show that this redshift-independent SFR/ ˙ M h -M h relation can reproduce the UV luminosity functions at z glyph[greaterorsimilar] 4 (e.g., Mason et al. 2015a; Harikane et al. 2018a; Tacchella et al. 2018; Bouwens et al. 2021) and the trend of the redshift evolution of the cosmic SFR density (e.g., Mason et al. 2015a; Harikane et al. 2018a; Tacchella et al. 2018; Oesch et al. 2018), a.k.a the cosmic star formation history or the Lilly-Madau plot (e.g., Lilly et al. 1996; Madau et al. 1996; Sawicki et al. 1997; Steidel et al. 1999; Bouwens et al. 2015; see review by Madau & Dickinson 2014). As discussed in Harikane et al. (2018a), this suggests a simple picture that the evolution of the cosmic SFR density is primarily driven by the steep increase of the number density of halos (and galaxies) due to the structure formation to z ∼ 4 -2, and the decrease of the accretion rate from z ∼ 2 to z ∼ 0 due to the cosmic expansion. However, the SFR/ ˙ M h -M h relation is only constrained at z ∼ 4 -7, and it is not known whether the relation evolves from z ∼ 4 to z ∼ 1 -3 or not, where the cosmic SFR density reaches its peak.  \nIn this work, we present new measurements of the restframe UV luminosity functions at z ∼ 4 -7 and clustering at z ∼ 2 -6 based on wide and deep optical images obtained in the Hyper-Suprime Cam Subaru Strategic Program (HSC-SSP) survey (Aihara et al. 2018, see also Miyazaki et al. 2012, 2018; Komiyama et al. 2018; Furusawa et al. 2018) and the CFHT large area U-band deep survey (CLAUDS; Sawicki et al. 2019). This paper is one in a series of papers from twin programs dedicated to scientific results on high redshift galaxies based on the HSCSSP survey data. One program is our luminous Lyman break galaxy (LBG) or dropout galaxy studies, named  \nTable 1 HSC-SSP Data Used for the z ∼ 4 -7 Selection  \nNote . -(1) Field name. (2) Right ascension. (3) Declination. (4) Effective area in deg 2 . (5)-(9) 5 σ limiting magnitude measured in 1 . '' 5 diameter circular apertures in the g -, r -, i -, z -, and y -bands. These limiting magnitudes are not corrected to total.  \nGreat Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH; Ono et al. 2018, Harikane et al. 2018a, Toshikawa et al. 2018). The other program is high redshift Ly α emitter (LAE) studies using HSC narrowband filters, Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH; Ouchi et al. 2018; Shibuya et al. 2018a,b; Konno et al. 2018; Harikane et al. 2018b, 2019; Inoue et al. 2018; Higuchi et al. 2019; Kakuma et al. 2019; Ono et al. 2021). Our new LBG catalogs are made public on our project webpage 1 or zenodo. 2  \nThis paper is organized as follows. We show the observational datasets in Section 2 and describe sample selections in Section 3. The results of the UV luminosity functions and clustering analysis are presented in Sections 4 and 5, respectively. We discuss our results in Section 6, and summarize our findings in Section 7. Throughout this paper we use the Planck cosmological parameter sets of the TT, TE, EE+lowP+lensing+ext result (Planck Collaboration et al. 2016): Ω m = 0 . 3089, Ω Λ = 0 . 6911, Ω b = 0 . 049, h = 0 . 6774, and σ 8 = 0 . 8159. We define r 200 that is the radius in which the mean enclosed density is 200 times higher than the mean cosmic density. To define the halo mass, we use M 200 that is the dark matter and baryon mass enclosed in r 200 . Note that this definition is the same as Harikane et al. (2016) but different from the one in Harikane et al. (2018a) who use the total dark matter mass without baryons. We assume the Salpeter (1955) initial mass function (IMF). All magnitudes are in the AB system (Oke & Gunn 1983).", '2.1. Subaru/HSC Data': "We use the internal S18A data release product taken in the HSC-SSP survey (Aihara et al. 2018) from 2014 March to 2018 January, which is basically identical to the version of the Public Data Release 2 (Aihara et al.  \n2019). 3 The HSC-SSP survey obtains deep optical imaging data with the five broadband filters, g , r , i , z , and y (Kawanomoto et al. 2018), which are useful to select z ∼ 4 -7 galaxies with the dropout selection technique. The HSC-SSP survey has three layers, the UltraDeep, Deep, and Wide, with different combinations of area and depth. Total effective survey areas of the data we use are ∼ 3, ∼ 18, and ∼ 288 deg 2 for the UltraDeep, Deep, and Wide layers, respectively (Table 1). Here we define the effective survey area as area where the number of visits in g , r , i , z , and y -bands are equal to or larger than threshold values after masking interpolated, saturated, or bad pixels, cosmic rays, and bright source halos (Coupon et al. 2018). The applied flags and threshold values are summarized in Table 2. In addition to these flags, we mask some regions that are affected by bright source halos or ghosts of bright sources but not flagged.  \nThe HSC data are reduced by the HSC-SSP collaboration with hscPipe (Bosch et al. 2018) that is the HSC data reduction pipeline based on the Large Synoptic Survey Telescope (LSST) pipeline (Ivezic et al. 2008; Axelrod et al. 2010). hscPipe performs all the standard procedures including bias subtraction, flat fielding with dome flats, stacking, astrometric and photometric calibrations, flagging, source detections and measurements, and construction of a multiband photometric catalog. The astrometric and photometric calibration are based on the data of Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) 1 imaging survey (Magnier et al. 2013; Schlafly et al. 2012; Tonry et al. 2012). PSFs are calculated in hscPipe , and typical full width at half maximum (FWHM) of the PSFs are 0 . '' 6 -0 . '' 9.  \nWe use forced photometry, which allows us to measure fluxes in multiple bands with a consistent aperture defined in a reference band. The reference band is i by default and is switched to z ( y ) for sources with no detection in the i ( z ) and bluer bands. In previous studies based on the S16A data release (e.g., Ono et al.  \nTable 2 Selection Criteria for Our Catalog Construction  \n2018; Harikane et al. 2018a; Toshikawa et al. 2018), the CModel magnitude (Abazajian et al. 2004) was used to measure total fluxes and colors of sources. However, we have found that some objects in the S18A data release have unnaturally bright CModel magnitudes compared to their aperture magnitudes, as also reported in Hayashi et al. (2020). Thus in this paper, we instead use magnitudes measured with a fixed 2 '' -diameter aperture after aperture correction, convolvedflux 0 20 , for measuring total fluxes and colors of sources. The aperture correction factor is calculated in each band assuming the point spread function. Among several magnitudes with different aperture sizes and corrections calculated with hscPipe , we have found that this magnitude provides the best match to the CModel magnitude in the S16A data release with the typical difference less than 5%. Limiting magnitudes and source detections are evaluated with magnitudes measured in a 1 . '' 5-diameter aperture, m aper . All the magnitudes are corrected for Galactic extinction (Schlegel et al. 1998). We measure the 5 σ limiting magnitudes which are defined as the 5 σ levels of sky noise in a 1 . '' 5-diameter aperture. The sky noise is calculated from fluxes in sky apertures which are randomly placed on the images in the reduction process. The limiting magnitudes measured in g , r , i , z , and y -bands are presented in Table 1.  \nWe select isolated or cleanly deblended sources from detected source catalogs available on the database that are provided by the HSC-SSP survey team. We require that none of the central 3 × 3 pixels are saturated, and there are no bad pixels in their footprint, like the definition of the effective area described above. We also require that there are no problems in measuring the CModel fluxes in the gri images for g -dropouts, in the riz images for r -dropouts, in the izy images for i -dropouts, and in the zy images for z -dropouts, except for the problem of unnaturally bright magnitudes described above. In addition, we remove sources if there are any problems in measuring their centroid positions in the ri images for g -dropouts, in the iz images for r -dropouts, in the zy im-  \nges for i -dropouts, and in the y image for z -dropouts. To remove severely blended sources, we apply a blendedness parameter threshold of b < 0 . 2 in the ri , iz , and zy -bands at z ∼ 4, 5, and 6, respectively These selection criteria are summarized in Table 2.", '2.2. CLAUDS Data': 'In the UltraDeep and Deep layers of the HSC-SSP survey, deep U -band images are taken in CLAUDS (Sawicki et al. 2019). These U -band images are useful to select z ∼ 2 -3 galaxies by using the U -dropout or BX/BM selection techniques (e.g., Steidel et al. 2003, 2004; Adelberger et al. 2004). The U -band images are obtained with two filters, u and u ∗ , because CFHT updated the MegaCam filter set by replacing the old u ∗ -filter by the new u -filter during the CLAUDS observing campaign (2014B to 2016B). Specifically, we have deep u ∗ -band images in the UD-SXDS, UD-COSMOS, and D-XMM-LSS fields, and deep u -band images in D-COSMOS, D-DEEP2-3 and D-ELAIS-N1 fields. Sawicki et al. (2019) describe the data reduction and procedures for making combined CLAUDS+HSC-SSP catalogs in detail. The 5 σ depths of the u ∗ and u -band images are typically 27.9 and 27.5 mag, respectively, sufficiently deep to select z ∼ 2 -3 galaxies.', '2.3. Spectroscopic Data': "We carried out spectroscopic follow-up observations for sources in our dropout catalogs at z ∼ 4 -7 with DEep Imaging Multi-Object Spectrograph (DEIMOS; Faber et al. 2003) on the Keck Telescope on 2018 August 11 (S18B-014, PI: Y. Ono), AAOmega+2dF (Sharp et al. 2006; Lewis et al. 2002) on the Anglo-Australian Telescope (AAT) from December 31 2018 to January 3 2019 (A/2018B/03, PI: Y. Ono), and the Faint Object Camera and Spectrograph (FOCAS: Kashikawa et al. 2002) on the Subaru Telescope on 2019 May 13 (S19A-006, PI: Y. Ono).  \nIn the DEIMOS observations, we used the 600ZD grating with the GG455 filter. The spectroscopic observa-  \ntions were made in multi-object slit mode. We used a total of two masks. Slit widths were 0. '' 8, and the integration time was 3600-6000 seconds per each mask. The DEIMOS spectra were reduced with the spec2d IDL pipeline developed by the DEEP2 Redshift Survey Team (Cooper et al. 2012; Newman et al. 2013). Wavelength calibrations was conducted by using the arc lamp emission lines. The spectral resolutions in an FWHM based on the widths of night-sky emission lines were ∼ 3 . 7 ˚ A. Flux calibration was achieved with data of a standard star G191B2B. The details of the DEIMOS observations will be presented in Ono et al. in prep.  \nIn the AAOmega+2dF observations, we used the X5700 dichroic beam splitter, the 580V grating with the central wavelength at 4821 ˚ A in the blue channel, and the 385R grating with the central wavelength at 7251 ˚ A in the red channel. This configuration covered a wavelength range of 3800 -8800 ˚ A with a resolution of R ∼ 1400. We used a total of four masks, two covering the UD-COSMOS field and two covering the UD-SXDS field. The integration time is 1800-7800 seconds per each mask, although weather conditions were not excellent. The spectra are reduced in the standard manner by using the OzDES pipeline (Yuan et al. 2015; Childress et al. 2017; Lidman et al. 2020).  \nIn the FOCAS observations, we used the VPH900 grism with the SO58 order-cut filter. The spectroscopic observations were made in the multi-object slit mode. We used a total of two masks. Slit widths were 0. '' 8, and the integration time was 7200 seconds per each mask. The FOCAS data were reduced with the focasred pipeline. Wavelength calibrations was conducted by using night-sky emission lines. The spectral resolution in an FWHM of FOCAS VPH900 based on the night-sky lines was ∼ 5 . 7 ˚ A. Flux calibration was performed with data of a standard star BD+28d4211. The details of the FOCAS observations will be presented in Ono et al. in prep.  \nA total of 55 dropout candidates were targeted in these observations. Target priorities were determined by apparent magnitudes, and brighter sources were assigned higher priorities. The apparent magnitude range of the targets was 19 -25 mag and most of them were 21 -24 mag.  \nIn addition to the observations described above, we include results of our observations with the Inamori Magellan Areal Camera and Spectrograph (IMACS; Dressler et al. 2011) on the Magellan I Baade telescope in 2007 - 2011 (PI: M. Ouchi). The IMACS observations were carried out on 2007 November 11-14, 2008 November 29-30, December 1-2, December 18-20, 2009 October 11-13, 2010 February 8-9, July 9-10, and 2011 January 3-4. In these observations, main targets were high redshift LAE candidates found in the deep Subaru SuprimeCam narrowband images obtained in the SXDS (Ouchi et al. 2008, 2010) and COSMOS fields (Murayama et al. 2007; Shioya et al. 2009). High redshift dropout galaxy candidates selected from deep broadband images in these two fields (Furusawa et al. 2008; Capak et al. 2007) were also observed as mask fillers. The data were reduced with the Carnegie Observatories System for MultiObject  \nSpectroscopy ( cosmos ) pipeline. 4", '3.1. Source Selection at z ∼ 4 -7': "First, to identify secure sources, we select sources whose signal-to-noise (SN) ratios are higher than 5 within 1 . '' 5-diameter apertures in the i -band for g -dropouts and in the z -band for r -dropouts. For i -dropouts, we select sources with signal-to-noise ratios higher than 5 and 4 in the z and y -bands, respectively, because the y -band images are relatively shallow. For the z -dropouts, we select sources with signal-to-noise ratios higher than 5 in the y -band. We then select dropout galaxy candidates by using their broadband spectral energy distribution (SED) colors. Like our previous studies (Ono et al. 2018; Harikane et al. 2018a; Toshikawa et al. 2018), we adopt the following color criteria:  \nFrom the source catalogs made in Section 2.1, we construct z ∼ 4 -7 dropout candidate catalogs based on the Lyman break color selection technique (e.g., Steidel et al. 1996; Giavalisco 2002). As shown in Figure 1, galaxy candidates can be selected based on their gri , riz , izy , and zy colors at z ∼ 4, 5, 6, and 7, respectively.  \ng -dropouts ( z ∼ 4)  \n( g -r > 1 . 0) ∧ (1)  \n( g -r > 1 . 5( r -i ) + 0 . 8) , (3)  \n( r -i < 1 . 0) ∧ (2)  \nr -dropouts ( z ∼ 5)  \n( r -i > 1 . 2) ∧ (4)  \n( r -i > 1 . 5( i -z ) + 1 . 0) , (6)  \n( i -z < 0 . 7) ∧ (5)  \ni -dropouts ( z ∼ 6)  \n( i -z > 1 . 5) ∧ (7)  \n( i -z > 2 . 0( z -y ) + 1 . 1) , (9)  \n( z -y < 0 . 5) ∧ (8)  \nz -dropouts ( z ∼ 7)  \nz -y > 1 . 6 . (10)  \nAs shown in Figure 1, these color criteria are set to avoid low redshift galaxies and stellar contaminants. Although only the z -y color is used in the z -dropout selection, we can efficiently select z ∼ 7 sources with this strict color criterion, as shown in previous studies similarly using the z -y color (e.g., Ouchi et al. 2009).  \nTo remove foreground interlopers, we exclude sources with continuum detections at > 2 σ levels in the g -band for r -dropouts, in the g - or r -bands for i -dropouts, and in the g , r , or i -bands for z -dropouts, using the 1 . '' 5 diameter aperture magnitudes. Since our z -dropout candidates are detected only in y -band images, we carefully check coadd and single epoch observation images of the selected candidates to remove spurious sources and moving objects.  \nUsing the selection criteria described above, we select a total of 1,978,462 dropout candidates at z ∼ 4 -7, consisting of 1,836,244 g -dropouts, 139,359 r -dropouts,  \n---Figure 1. Two-color diagrams to select dropout sources. The left, middle, and right panels show two-color diagrams to select g -dropout ( z ∼ 4), r -dropout ( z ∼ 5), and i -dropout ( z ∼ 6) sources, respectively. The red lines indicate color criteria we use to select dropout sources (Equations (1)-(9)), and the red circles are spectroscopically-identified sources in the UltraDeep layers (for the left panel) and in all layers (for the middle and right panels). The solid black lines are colors of star forming galaxies expected from the Bruzual & Charlot (2003) model as a function of redshift. As model parameters, we adopt the Salpeter (1955) IMF, an age of 70 Myr after the initial star formation, metallicity of Z/Z glyph[circledot] = 0 . 2, and the Calzetti et al. (2000) dust extinction law with reddening of E ( B -V ) = 0 . 16 (see Section 3.2). The circles on the line show their redshifts with an interval of ∆ z = 0 . 1. The blue circles are z = 0 -3 sources spectroscopically identified in the UD-COSMOS region. The dotted, dashed, and dot-dashed lines are, respectively, typical spectra of elliptical, Sbc, and irregular galaxies (Coleman et al. 1980) redshifted from z = 0 to z = 2 (for the left panel) and to z = 3 (for the middle and right panels). The black stars indicate Galactic stars taken from Gunn & Stryker (1983) and L /T dwarfs from Knapp et al. (2004).  \n<!-- image -->  \nTable 3 Number of Sources in Galaxy Samples Used in This Work.  \n- † Selected in C. Liu et al. in prep.  \n∗ Selected in This work.  \n2,567 i -dropouts, and 292 z -dropouts. Our sample is selected from the 307 . 9 deg 2 wide area data corresponding to a 5 . 92 Gpc 3 survey volume, and is the largest sample of the high-redshift ( z glyph[greaterorsimilar] 4) galaxy population to date. Especially, combined with z ∼ 2 -3 galaxies selected later in Section 3.5, we have a total of 4 , 100 , 221 galaxies at z ∼ 2 -7, which is the largest among high redshift galaxy studies. Table 3 summarizes the number of dropout candidates in each field, and Figure 2 shows examples of sky distributions of the dropouts. The differences in the numbers of the selected candidates mainly come from the differences in the survey areas and depths.", '3.2. Spectroscopically-Identified Sources': "In our samples, a total of 46 sources are identified as z > 3 objects through our spectroscopic  \nfollow-up observations with DEIMOS, AAOmega+2dF, and FOCAS (Section 2.3). Redshifts are determined based on the Ly α line and/or Lyman break. Figure 3 shows examples of spectroscopically identified sources, HSC J160953+532821, HSC J161207+555919, HSC J020834 -021239, and HSC J022552 -054439. HSC J160953+532821 is a faint quasar at z = 6 . 923 with a UV magnitude of M UV = -22 . 7 mag, confirmed in our FOCAS observations. This source is also identified in Matsuoka et al. (2019). An FWHM of its Ly α line is ∼ 2500 km s -1 , and a rest-frame Ly α equivalent width (EW) is EW 0 Ly α ∼ 600 ˚ A. HSC J161207+555919 is a bright ( M UV = -23 . 0 mag) galaxy at z = 6 . 788 with a narrow Ly α emission line of < 300 km s -1 and EW 0 Ly α ∼ 3 ˚ A, confirmed in our FOCAS observations. This source  \nFigure 2. Examples of sky distributions of dropout sources. The red squares, orange circles, and black dots are, respectively, the positions of sources whose magnitudes are i < 22 . 0, 22 . 0 -22 . 5, and 22 . 5 -25 . 0 (24 . 0) mag for z ∼ 4 sources, z < 22 . 5, 22 . 5 -23 . 0, and 23 . 0 -25 . 5 (24 . 5) mag for z ∼ 5 sources, and y < 23 . 0, 23 . 0 -23 . 5, and 23 . 5 -25 . 5 (25 . 0) mag for z ∼ 6 and 7 sources, in the UltraDeep (Deep and Wide) layers. The scale on the map is marked in degrees and in the projected distance (comoving Gpc).  \n<!-- image -->  \nFigure 3. Spectra of HSC J160953+532821 at z = 6 . 923, HSC J161207+555919 at z = 6 . 788, HSC J020834 -021239 at z = 4 . 088, and HSC J022552 -054439 at z = 3 . 647 from top to bottom, obtained in our spectroscopic follow-up observations. In each figure, the top panel shows the two-dimensional spectrum (black is positive) and the bottom panel shows the one-dimensional spectrum. In the top panel, our dropout galaxy is located at the center in the spatial direction, and the spatial range is ± 3 '' . In the bottom panel, the black line indicates the spectrum of the object, and the blue line shows the sky spectrum with an arbitrary normalization. For the sources with weak Ly α emission, we also plot the averaged spectra over 200 ˚ A bins with red filled circles to show the Lyman break features.  \n<!-- image -->  \nWavelength  \nλ  \nobs  \n[  \nÅ  \n]  \nis identified as a narrow-line quasar in Matsuoka et al. (2019). The Lyman break feature can be seen in the spectrum, whose redshift is consistent with that derived from the Ly α emission line. HSC J020834 -021239 and HSC J022552 -054439 are bright galaxies at z = 4 . 088 and 3 . 647 identified in our DEIMOS observations, with M UV = -22 . 7 and -21 . 8 mag, FWHMs of 250 and 180 km s -1 , and EW 0 Ly α ∼ 40 and 9 ˚ A, respectively.  \nIn addition, we incorporate results of our spectroscopic observations for high redshift galaxies with Magellan/IMACS. We also check spectroscopic catalogs in other studies (Cuby et al. 2003, Ouchi et al. 2008, Saito et al. 2008, Willott et al. 2009, Willott et al. 2010, CurtisLake et al. 2012, Mallery et al. 2012, Masters et al. 2012, Le F'evre et al. 2013, Willott et al. 2013, Kriek et al. 2015, Ba˜nados et al. 2016, Hu et al. 2016, Matsuoka et al. 2016, Momcheva et al. 2016, Toshikawa et al. 2016, Wang et al. 2016, Hu et al. 2017, Jiang et al. 2017, Masters et al. 2017, Tasca et al. 2017, Yang et al. 2017, Hasinger et al. 2018, Matsuoka et al. 2018a, Matsuoka et al. 2018b, Ono et al. 2018, Pˆaris et al. 2018, Pentericci et al. 2018, Shibuya et al. 2018b, Matsuoka et al. 2019, Harikane et al. 2020b, Harikane et al. 2020a, Zhang et al. 2020, Endsley et al. 2021, and Garilli et al. 2021). We adopt their classifications between galaxies and AGNs in their catalogs, which are mostly based on apparent AGN features such as broad emission lines. For the catalogs of the VIMOS VLT Deep Survey (VVDS; Le F'evre et al. 2013) and the VIMOS Ultra Deep Survey (VUDS; Tasca et al. 2017), we take into account sources whose reliabilities of the redshift determinations are > 70 -75%, i.e., sources with redshift reliability flags of 2, 3, 4, 9, 12, 13, 14, and 19. Here we focus on sources with spectroscopic redshifts of z spec > 3 in these catalogs.  \nIn total, 1037 dropouts in our sample have been spectroscopically identified at z spec ≥ 3 in our observations and previous studies, including 770 galaxies and 267 AGNs. These sources are listed in Table 10, and the redshift distributions of the sources are shown in Figure 4. Figure 1 shows the distributions of the spectroscopically identified sources at z spec > 3 in the two-color diagrams. We also plot sources in the UD-COSMOS field with spectroscopic redshifts of z spec < 3 as foreground interlopers. In addition, the tracks of model spectra of young star-forming galaxies that are produced with the stellar population synthesis code galaxev (Bruzual & Charlot 2003) are shown. As model parameters, the Salpeter (1955) IMF, an age of 70 Myr after the initial star formation, and metallicity of Z/Z glyph[circledot] = 0 . 2 are adopted. We use the Calzetti et al. (2000) dust extinction law with reddening of E ( B -V ) = 0 . 16. The IGM absorption is considered following the prescription of Madau (1995). The colors of the spectroscopically identified galaxies are broadly consistent with those expected from the model spectra.", '3.3. Selection Completeness and Redshift Distribution': 'The selection completeness and redshift distributions of dropout candidates at z ∼ 4 -7 are estimated based on results of Monte Carlo simulations in Ono et al. (2018). Ono et al. (2018) run a suite of Monte Carlo simulations with an input mock catalog of high redshift galaxies with the size distribution of Shibuya et al. (2015), the Sersic  \nFigure 4 shows results of the selection completeness estimates as a function of redshift. The average redshift values are roughly ¯ z = 3 . 8 for g -dropouts, ¯ z = 4 . 9 for r -dropouts, ¯ z = 5 . 9 for i -dropouts, and ¯ z = 6 . 9 for z -dropouts. In Figure 4, we also show the redshift distributions of the spectroscopically identified galaxies in our samples (Section 3.2). The redshift distributions of the spectroscopically identified sources are broadly consistent with the results of our selection completeness simulations. However, the distributions of the spectroscopically identified sources in the g -, r -, and i -dropout samples appear to be shifted toward slightly higher redshifts compared to the simulation results. This is probably because the spectroscopically identified sources are biased to ones with strong Ly α emission in the g - and r -dropout samples, and ones identified in the SHELLQs project searching for z ∼ 6 -7 quasars (e.g., Matsuoka et al. 2016) in the i -dropout sample. In particular, the redshift distribution of the spectroscopically identified r -dropouts has a secondary peak at z ∼ 5 . 7, which is  \n<!-- image -->  \nRedshift  \nFigure 4. Selection completeness estimates for our z ∼ 4, 5, 6, and 7 samples. The black curves correspond to the results of the Monte Carlo calculations in Section 3.3 averaged over the Wide, Deep, and UltraDeep layers. Average redshifts of these samples are 3.8, 4.9, 5.9, and 6.9. The histograms show redshift distributions of spectroscopically identified sources in the z ∼ 4 (blue), z ∼ 5 (green), z ∼ 6 (orange), and z ∼ 7 (red) samples.  \nindex of n = 1 . 5, the intrinsic ellipticities of 0 . 0-0 . 8. To produce galaxy SEDs, the stellar population synthesis model of galaxev (Bruzual & Charlot 2003) is used with the Salpeter (1955) IMF, a constant rate of star formation, age of 25 Myr, metallicity of Z/Z glyph[circledot] = 0 . 2, and the Calzetti et al. (2000) dust extinction ranging from E ( B -V ) = 0 . 0-0 . 4, corresponding to the UV spectral slope of -3 glyph[lessorsimilar] β UV glyph[lessorsimilar] -1. We do not include Ly α emission in the galaxy SEDs because the line fluxes are typically not significant compared to the continuum in the broadband fluxes. The IGM attenuation is taken into account by using the prescription of Madau (1995). Different simulations are carried out for the Wide, Deep, and UltraDeep layers by using the SynPipe software (Huang et al. 2018), which utilizes GalSim v1.4 (Rowe et al. 2015) and the HSC pipeline hscpipe (Bosch et al. 2018). We insert large numbers of artificial sources into HSC images. Then we select high redshift galaxy candidates with the same selection criteria, and calculate the selection completeness as a function of magnitude and redshift, C ( m,z ), averaged over UV spectral slope β UV weighted with the β UV distribution of Bouwens et al. (2014). Then the obtained completeness is scaled based on the limiting magnitudes to correct for differences in depths of the S18A data in this study and S16A data used in Ono et al. (2018).  \ncaused by z = 5 . 7 Ly α emitters found in the Subaru Suprime-Cam and HSC narrowband surveys in the literature (e.g., Ouchi et al. 2008; Shibuya et al. 2018b). Another possible reason is the systematic uncertainty of the IGM attenuation model in the simulations. The model of Madau (1995) predicts lower redshift dropouts than Inoue et al. (2014), which may explain the discrepancy. For the z -dropout sample, the redshift distribution of the spectroscopically identified sources is shifted to the lower redshift, because of the increasing fraction of the neutral hydrogen at z > 7, which resonantly scatters Ly α photons.', '3.4. Contamination': "Some foreground objects such as red galaxies at low redshifts can satisfy our color selection criteria due to photometric scatters, although intrinsically they do not enter the color selection window. This happens especially in the Wide and Deep layers, whose limiting magnitudes are relatively shallow. We evaluate contamination fractions in our dropout samples with the following three methods.  \nThe first method is the one using estimates of photometric redshifts. We use photometric redshifts estimated with the MIZUKI code (Tanaka 2015; Tanaka et al. 2018). Here a foreground interloper is defined as a source whose 95% upper bound of the photometric redshift is less than z = 2. We derive the fraction of the foreground interlopers as a function of the i - ( z -) band magnitude in our z ∼ 4 ( z ∼ 5) dropout samples in the Wide and Deep layers. The derived contamination fractions are presented in Table 4. In the UltraDeep layer, the fractions of the interlopers in our z ∼ 4 -5 samples are negligibly small ( < 10%) at m UV ≥ 25 mag. The contamination fraction becomes higher for brighter sources at glyph[lessorsimilar] 24 mag, because the number density of high redshift galaxies decreases to brighter magnitudes, while that of foreground interlopers (e.g., low redshift passive galaxies) does not significantly change (e.g., Muzzin et al. 2013). We do not derive the contamination fraction of the z ∼ 6 -7 sources with the MIZUKI code, because accuracies of the photometric redshifts are not high due to the limited number of available bands redder than the Lyman break. In the UD-COSMOS field, we also check photometric redshifts in the COSMOS 2015 catalog (Laigle et al. 2016) that are determined with multi-band photometric data including near-infrared images such as Spizter /IRAC, useful to eliminate stellar contaminants (Stevans et al. 2018). We match our dropout sources with the COSMOS 2015 catalog within a 1 '' search radius of the object coordinate. For bright z ∼ 4 -5 sources detected with > 10 σ significance levels in the HSC i -band, whose flux measurements are reliable, typically less than 20% of them are classified as foreground interlopers, consistent with the estimates of the MIZUKI code. For fainter sources at z ∼ 4 -5, about ∼ 70% of them are classified as z > 3 galaxies. The contamination rate for z ∼ 6 -7 sources is less than 40%, although the number of sources is small and their flux measurements have relatively large uncertainties due to their faintness compared to the bright z ∼ 4 -5 sources. Stellar templates are also fitted for the sources in the COSMOS 2015 catalog, and only 3.5% of them are classified as stars with χ 2 star < χ 2 galaxy (see  \nTable 4 Estimated Contamination Fraction for the z ∼ 4 and z ∼ 5 Samples in UltraDeep (UD), Deep (D), and Wide (W) Layers  \nNote . - These contamination fractions are estimated based on the photometric redshift analysis (the first method in Section 3.4).  \nLaigle et al. 2016). We also check the COSMOS Hubble ACS catalog (Leauthaud et al. 2007), and only 3.7% of them are point sources with the SExtractor (Bertin & Arnouts 1996) stellarity parameters of > 0 . 9 (1=star; 0=galaxy). These analyses indicate that the fraction of the stellar contamination is negligibly small.  \nThe second method is the one using the spectroscopic redshift catalog created in Section 2.3. We estimate the contamination fraction in the z ∼ 4 sample using the results in our AAOmega+2dF spectroscopy and the VVDS spectroscopy, which target a sufficient number of bright z ∼ 4 sources, and whose target selections are not significantly biased to low or high redshift sources. Foreground interlopers are identified based on continuum emission bluewards of the expected wavelength of the Lyman break, and/or rest-frame optical emission lines. At z glyph[greaterorsimilar] 5, we cannot derive robust contamination fractions because of the small number of spectroscopically confirmed sources in the AAOmega+2dF and VVDS data. Nonetheless, we find that a total of 80 sources from our z ∼ 6 -7 dropout samples ( y = 21 . 0 -25 . 6 mag) are spectroscopically identified in the entire spectroscopic catalog, and all of them are at z spec > 5, although it is possible that the actual contamination rate is higher than inferred from these numbers due to various biases including the publication bias and the fact that spectroscopic observations usually prioritize the most promising candidates.  \nThe third method is a simulation with shallower data, in the same manner as Ono et al. (2018). We use a shallower dataset whose depth is comparable with that of the Wide layer in the UD-COSMOS field, the Widelayer-depth COSMOS data. We assume that the UDCOSMOS data are sufficiently deep and the contamination fraction in our dropout selections is small. First, we select objects that do not satisfy our selection criteria at each redshift from the UD-COSMOS catalog. Then, we find the closest source in the Wide-layer-depth COSMOS catalog that matches within a 1 '' search ra-  \nof the object coordinate. If the objects satisfy our selection criteria for the Wide-layer dropout, we regard them as foreground interlopers, and calculate their number densities. Based on comparisons between the surface number densities of interlopers and those of the selected dropouts, we estimate the fractions of foreground interlopers that satisfy our color selection criteria due to the photometric scatters. The estimated contamination fractions are ∼ 0 -40% for sources with 24 -25 mag at z ∼ 4 -5, comparable to those in Ono et al. (2018). For the z ∼ 6 -7 dropout samples, we cannot estimate the surface number densities of interlopers by adopting this method, because the number densities of such sources in the shallower depth COSMOS field data are too small due to the limited survey area of the UD-COSMOS field.  \nWe find that the three methods above give the contamination fractions consistent with each other within their uncertainties at z ∼ 4 -5. As the contamination fractions used in derivations of luminosity functions later, we adopt the fractions determined based on the photometric redshifts, given their high accuracies compared to those of the other methods. For the z ∼ 4 -5 samples in the UltraDeep layer and the z ∼ 6 -7 samples, we assume that contamination fractions are negligibly small based on the results of the photometric and spectroscopic redshifts.", '3.5. Source Selection at z ∼ 2 -3': 'In addition to the z ∼ 4 -7 catalogs, we also use catalogs of galaxy candidates at z ∼ 2 -3 to study clustering properties. We use BM, BX, and U -dropout galaxy catalogs at z ∼ 1 . 7, 2 . 2, and 3 constructed in C. Liu et al. in prep. Here we briefly describe our source selection at z ∼ 2 -3. The z ∼ 2 -3 galaxy candidates are selected from a combined CLAUDS+HSC-SSP catalog made by Sawicki et al. (2019) based on hscpipe . Note that the HSC-SSP data in the combined CLAUDS+HSCSSP catalog is based on the S16A internal data release, which is different from the S18A data release product we use for the z ∼ 4 -7 selection. As described in Section 2.2, the U-band images are obtained with two filters, the u and u ∗ -bands. For the u ∗ -band filter, we adopt color criteria same as Hildebrandt et al. (2009), who select z ∼ 3 dropout galaxies with the similar filter set to this study:  \nU -dropouts ( z ∼ 3)  \n( g -r < 1 . 2) ∧ (11)  \n( u ∗ -g > 1 . 5( g -r ) + 0 . 75) . (13)  \n( u ∗ -g > 0 . 9) ∧ (12)  \nFor the BX and BM galaxies, we adopt the following color criteria, respectively: BX ( z ∼ 2 . 2)  \n( g -r > -0 . 5) ∧ (14)  \n( u ∗ -g > 0 . 42( g -r ) + 0 . 54) ∧ (16)  \n( u ∗ -g > 2 . 57( g -r ) + 0 . 21) ∧ (15)  \n(( u ∗ -g < 0 . 9) ∨ ( u ∗ -g < 1 . 5( g -r ) + 0 . 75)) , (17)  \nBM ( z ∼ 1 . 7)  \n( u ∗ -g > -2 . 5( g -r ) + 0 . 98) ∧ (19)  \n( g -r > -0 . 5) ∧ (18)  \n( u ∗ -g > 0 . 93( g -r ) + 0 . 12) ∧ (20)  \n( u ∗ -g < 0 . 42( g -r ) + 0 . 54) . (21)  \nFor the ugr filter set, we define selection criteria by comparing the positions of stars in the ugr and u ∗ gr diagrams:  \nU -dropouts ( z ∼ 3)  \n( g -r < 1 . 2) ∧ (22)  \n( u -g > 1 . 99 × g -r +0 . 68) , (24)  \n( u -g > 0 . 98) ∧ (23)  \nBX ( z ∼ 2 . 2)  \n( g -r > -0 . 5) ∧ (25)  \n( u -g > 0 . 42( g -r ) + 0 . 54) ∧ (27)  \n( u -g > 3 . 07( g -r ) + 0 . 14) ∧ (26)  \n(( u -g < 0 . 9) ∨ ( u -g < 1 . 5( g -r ) + 0 . 75)) . (28)  \nSelection criteria of BM galaxies for the ugr filter set are the same as those for the u ∗ gr filter set. Details of the selection are presented in C. Liu et al. in prep. A total of 935,804, 405,469, and 780,486 galaxy candidates are selected at z ∼ 1 . 7, 2 . 2, and 3, respectively. The number densities of the selected galaxies are comparable to previous studies. The selection completeness and contamination fraction of these samples are described in C. Liu et al. in prep.', '4. UV Luminosity Function': '4.1. Dropout UV Luminosity Function  \n4.1.1. Derivation  \nWe derive the rest-frame UV luminosity functions of z ∼ 4 -7 dropout sources by applying the effective volume method (Steidel et al. 1999). Based on the results of the selection completeness simulations in Section 3.3, we estimate the effective survey volume per unit area as a function of the apparent magnitude,  \nV eff ( m ) = ∫ C ( m,z ) dV ( z ) dz dz, (29)  \nwhere C ( m,z ) is the selection completeness, i.e., the probability that a galaxy with an apparent magnitude m at redshift z is detected and satisfies the selection criteria, and dV ( z ) /dz is the differential comoving volume as a function of redshift.  \nThe space number densities of the dropouts that are corrected for incompleteness and contamination effects are obtained by calculating  \nψ ( m ) = [1 -f cont ( m )] n raw ( m ) V eff ( m ) , (30)  \nwhere n raw ( m ) is the surface number density of selected dropouts in an apparent magnitude bin of m , and f cont ( m ) is the contamination fraction in the magnitude bin estimated in Section 3.4. The 1 σ uncertainties are calculated by taking account of the Poisson confidence limits (Gehrels 1986) on the numbers of the sources. To  \ncalculate 1 σ uncertainties of the space number densities of dropouts, we consider uncertainties of the surface number densities and the contamination fractions.  \nWe convert the number densities of dropouts as a function of apparent magnitude, ψ ( m ), into the UV luminosity functions, Φ[ M UV ( m )], which is the number densities of dropouts as a function of rest-frame UV absolute magnitude. We calculate the absolute UV magnitudes of dropout samples from their apparent magnitudes using their averaged redshifts ¯ z :  \nM UV = m +2 . 5 log(1 + ¯ z ) -5 log ( d L (¯ z ) 10 pc ) +( m UV -m ) , (31)  \nwhere d L is the luminosity distance in units of parsecs and ( m UV -m ) is the K -correction term between the magnitude at rest-frame UV and the magnitude in the bandpass that we use. We define the UV magnitude, m UV , as the magnitude at the rest-frame 1500 ˚ A. For the apparent magnitude m , we use a magnitude in a band whose central wavelength is the nearest to the restframe wavelength of 1500 ˚ A, namely i -, z -, y -, and y -bands for g -, r -, i -, and z -dropouts, respectively. We set the K -correction term to be 0 by assuming that dropout galaxies have flat UV continua, i.e., constant f ν in the rest-frame UV ( m UV = m ). Note that this assumption does not have a significant impact on the calculated UV magnitudes. If we vary the UV slope ( β UV ) with -2 . 5 < β UV < -1 . 5 (Bouwens et al. 2014), the calculated UV magnitude differs only within 0.1 mag.', '4.1.2. Results': 'The top panel of Figure 5 shows our derived luminosity function for g -dropouts at z ∼ 4 with previous studies. Our measurements have smaller error bars compared to Ono et al. (2018) because of the improved constraint on the contamination fraction (Section 3.4). Our measurements agree well with previous studies of quasars at M UV glyph[lessorsimilar] -24 (e.g., Akiyama et al. 2018), studies of galaxies at M UV glyph[greaterorsimilar] -22 (e.g., Bouwens et al. 2021; Finkelstein et al. 2015b), and studies of galaxies and AGNs (Ono et al. 2018; Stevans et al. 2018; Adams et al. 2020).  \nOur derived luminosity functions at z ∼ 5 and 6 are shown in the bottom panel of Figure 5 and the top panel of Figure 6, respectively. Similar to the z ∼ 4 result, our measurements agree well with previous studies of quasars at M UV glyph[lessorsimilar] -24 mag (e.g., Niida et al. 2020; Matsuoka et al. 2018c), studies of galaxies at M UV glyph[greaterorsimilar] -22 mag (e.g., Bouwens et al. 2021; Finkelstein et al. 2015b), and studies of galaxies and AGNs (Ono et al. 2018). At z ∼ 7, our derived luminosity function agrees with previous studies (e.g., Ono et al. 2018; Bowler et al. 2017), as shown in the bottom panel of Figure 6. Table 5 summarizes our measurements of the luminosity functions at z ∼ 4 -7.  \nThese agreements clearly indicate that the dropout luminosity function is a superposition of the AGN luminosity function (dominant at M UV < -24 mag) and the galaxy luminosity function (dominant at M UV > -22 mag). In our dropout selection, we probe redshifted Lyman break features of high redshift galaxies. However, high redshift AGNs also have similar Lyman break features. Thus it is expected that our dropout sample is composed of both galaxies and AGNs. Indeed, as described in Section 3.2, our dropout samples include both  \nspectroscopically identified galaxies and AGNs. Based on our spectroscopic results and the literature, we derive the galaxy fraction which is the number of spectroscopically confirmed high redshift galaxies divided by the sum of the numbers of spectroscopically confirmed galaxies and AGNs. The derived galaxy fractions for our z ∼ 4 -7 samples in each magnitude bin are presented in Figures 5 and 6. As shown in Figures 5 and 6, the galaxy fractions at z ∼ 4 -7 are about 0% at M UV < -24 mag, but then increase with increasing magnitude and reach about 100% at M UV > -22 mag. These results further suggest that our luminosity functions are dominated by AGNs at the bright end, and by galaxies at the faint end. The very wide area and deep depth of the HSC-SSP survey allow us to bridge the UV luminosity functions of high redshift galaxies and AGNs, both of which can be selected with redshifted Lyman break features (Ono et al. 2018).', '4.1.3. Fitting the Dropout Luminosity Functions': "We investigate the shape of the UV luminosity functions of dropout sources (galaxies+AGNs) by fitting them with several functional forms. Figure 7 shows our UV luminosity functions at z ∼ 4 -7 with previous results of galaxies based on the Hubble data (Bouwens et al. 2021; Ishigaki et al. 2018) and of quasars based on the HSC-SSP and Sloan Digital Sky Survey (SDSS) data (Akiyama et al. 2018; Niida et al. 2020; Matsuoka et al. 2018c). The combination of our results with the previous work reveals the UV luminosity functions in a very wide magnitude range of -29 glyph[lessorsimilar] M UV glyph[lessorsimilar] -14 mag, corresponding to the luminosity range of 0 . 002 L ∗ UV glyph[lessorsimilar] L UV glyph[lessorsimilar] 2000 L ∗ UV . We also show luminosity functions at z ∼ 3 taken from Moutard et al. (2020), Bouwens et al. (2021), and Zhang et al. (2021). As discussed in Section 4.1.2, the dropout luminosity function is a superposition of the AGN luminosity function and the galaxy luminosity function. Thus we simultaneously fit the AGN and galaxy luminosity functions. For the AGN luminosity function, we fit with a double-power law (DPL) function that is widely used in studies of AGNs:  \nφ ( L ) dL = φ ∗ [ ( L L ∗ ) -α + ( L L ∗ ) -β ] -1 dL L ∗ , (32)  \nwhere φ ∗ is the overall normalization, L ∗ is the characteristic luminosity, and α and β are the faint and brightend power-law slopes, respectively. We define a DPL function as a function of absolute magnitude, Φ( M UV ), as φ ( L ) dL = Φ( M UV ) dM UV ,  \nΦ( M UV ) = ln 10 2 . 5 φ ∗ × [ 10 0 . 4( α +1)( M UV -M ∗ UV ) +10 0 . 4( β +1)( M UV -M ∗ UV ) ] -1 , (33)  \nwhere M ∗ UV is the characteristic magnitude. For the galaxy luminosity function, we fit with a DPL function or the Schechter function (Schechter 1976):  \nφ ( L ) dL = φ ∗ ( L L ∗ ) α exp ( -L L ∗ ) d ( L L ∗ ) , (34)  \nwhere φ ∗ , L ∗ UV , and α are the overall normalization, the characteristic luminosity, and the faint power-law  \nm  \nUV  \nFigure 5. Each bottom panel shows rest-frame UV luminosity functions of dropout sources (including galaxies and AGNs) at z ∼ 4 and 5. The red circles show our results based on the HSC-SSP survey data. For comparison, we also show previous results for galaxies+AGNs in Ono et al. (2018, cyan circles), Stevans et al. (2018, orange squares), and Adams et al. (2020, orange diamonds), galaxies in Bouwens et al. (2021, magenta squares), Finkelstein et al. (2015b, green circles), van der Burg et al. (2010, blue triangles), and Pell'o et al. (2018, orange diamonds), and AGNs in Akiyama et al. (2018, black squares), Parsa et al. (2018, black triangles), Giallongo et al. (2015, black crosses), Glikman et al. (2011, black pluses), Niida et al. (2020, black diamonds), and Kim et al. (2020, black stars). Each top panel shows a fraction of galaxies in our dropout sample based on spectroscopic results. For the denominator of the fraction, the sum of the numbers of galaxies and AGNs is used. Note that this fraction is estimated based on various spectroscopic catalogs, and its uncertainty is discussed in Section 4.2.1.  \n<!-- image -->  \nm  \nUV  \nFigure 6. Same as Figure 5 but at z ∼ 6 and 7. For comparison, we also show previous results for galaxies+AGNs in Ono et al. (2018, cyan circles), galaxies in Bouwens et al. (2021, magenta squares), Pell'o et al. (2018, orange diamonds), Ishigaki et al. (2018, orange diamonds), Bowler et al. (2017, 2015, blue triangles), Finkelstein et al. (2015b, green circles), Schenker et al. (2013, brown triangles), and McLure et al. (2013, purple circles), and AGNs in Matsuoka et al. (2018c, black squares), Parsa et al. (2018, black triangles), Onoue et al. (2017, black stars), Jiang et al. (2016, black circles), Giallongo et al. (2015, black crosses), Kashikawa et al. (2015, black triangles), and Willott et al. (2010, black pluses). In the z ∼ 7 panel, the z ∼ 6 AGN luminosity functions are plotted because there are no results of the AGN luminosity function at z ∼ 7.  \n<!-- image -->  \nUV  \nTable 5 Obtained Dropout (Galaxy+AGN) and Galaxy UV luminosity functions at z ∼ 4, 5, 6, and 7.  \n† The galaxy luminosity function derived by using Equation (36) with the spectroscopic galaxy fractions presented in Figures 5 and 6.  \nslope, respectively. We define the Schechter function as a function of absolute magnitude, Φ( M UV ), as φ ( L ) dL = Φ( M UV ) dM UV ,  \nΦ( M UV ) = ln 10 2 . 5 φ ∗ 10 -0 . 4( M UV -M ∗ UV )( α +1) × exp ( -10 -0 . 4( M UV -M ∗ UV ) ) . (35)  \nFigure 7 shows the best-fit results in cases of the DPL+DPL and DPL+Schechter functions at z ∼ 3 -7. Note that in the fitting at z ∼ 7, we fix the parameters of the AGN luminosity function to values in Matsuoka et al. (2018c) with decreasing the φ ∗ parame-  \n0.7 dex following φ ∗ ∝ 10 -0 . 7( z -6) as assumed in Matsuoka et al. (2018c), because there are no measurements of the AGN luminosity function at z ∼ 7. The luminosity functions in the very wide luminosity range of -29 glyph[lessorsimilar] M UV glyph[lessorsimilar] -14 mag are well fitted with either the DPL+DPL or DPL+Schechter functions, as partly shown at z ∼ 4 in previous studies (Stevans et al. 2018; Adams et al. 2020), except for z ∼ 7 where the DPL+DPL functions provide a better fit. Table 6 summarizes the best-fit parameters and reduced χ 2 values of the two cases. We find that the best-fit parameters are roughly comparable to those of galaxy and AGN luminosity functions in previous studies (e.g., Ono et al. 2018;  \nFigure 7. Rest-frame UV luminosity functions of all rest-UV selected sources (galaxies+AGNs) at z ∼ 3, 4, 5, 6, and 7 from top to bottom. The red circles show our results based on the HSC-SSP survey data at z ∼ 4 -7 and results in Moutard et al. (2020) at z ∼ 3. The magenta squares and orange diamonds are results for galaxies taken from Bouwens et al. (2021) and Ishigaki et al. (2018), respectively. The orange squares are results for AGNs in Zhang et al. (2021), Akiyama et al. (2018), Niida et al. (2020), and Matsuoka et al. (2018c) at M UV < -24 mag. The blue (green) lines show the best-fit DPL+DPL (DPL+Schechter) functions. Note that at z ∼ 7, we fix the parameters of the AGN luminosity function to values in Matsuoka et al. (2018c) with decreasing the φ ∗ parameter by 0.7 dex following φ ∗ ∝ 10 -0 . 7( z -6) as assumed in Matsuoka et al. (2018c).  \n<!-- image -->  \nAkiyama et al. 2018; Niida et al. 2020; Matsuoka et al. 2018c).", '4.1.4. Redshift Evolution': 'Figure 8 summarizes UV luminosity function estimates at z ∼ 4 -7 in this work and the literature, and the bestfit DPL+DPL functions. We also plot the rest-frame UV luminosity functions at z ∼ 0 -3 from Moutard et al. (2020). Like our z ∼ 4 -7 results, the luminosity functions at z ∼ 0 -3 also show number density excesses at the bright end compared to the Schechter functions, which are dominated by AGNs (see discussions in Moutard et al. 2020). Indeed, the number densities of bright sources ( M UV glyph[lessorsimilar] -23 mag) at z ∼ 3 are comparable to the rest-frame UV luminosity function of spectroscopically identified AGNs in Zhang et al. (2021). Interestingly, the number density of typical galaxies ( M UV = M ∗ UV ∼ -21 mag) increases only by a factor of ∼ 3 from z ∼ 6 to 3, while the number density of typical quasars ( M UV = M ∗ UV ∼ -27 mag) significantly increases by a factor of ∼ 100, consistent with previous studies of the quasar luminosity functions (e.g., Matsuoka et al. 2018c; Niida et al. 2020). This indicates a very rapid growth of AGNs in the first 1.5 Gyr. If we extrapolate this evolution to higher redshift by assuming Φ ∝ 10 -0 . 7(1+ z ) (Matsuoka et al. 2018c), the number density of bright quasars ( M UV glyph[lessorsimilar] -26 mag) will be very small ( glyph[lessorsimilar] 2 × 10 -10 mag -1 Mpc -3 ) at z ≥ 7. More specifically, the number density of typical quasars increases by a factor of 10 from z ∼ 6 to 5, but increases only by a factor of 3 from z ∼ 4 to 3, indicating the accelerated evolution of the quasar luminosity function at z ∼ 3 -6 (Niida et al. 2020).', '4.2. Galaxy UV Luminosity Function 4.2.1. Derivation and Results': "We estimate the galaxy UV luminosity functions in a wide magnitude range by considering the contributions from AGNs in our dropout luminosity function measurements. To subtract the AGN contributions, we use the galaxy fraction estimates based on the spectroscopy shown in Figures 5 and 6. We multiply the dropout (galaxy+AGN) UV luminosity functions by the spectroscopic galaxy fractions, f galaxy , and obtain the galaxy luminosity functions, Φ galaxy :  \nΦ galaxy ( M UV ) = f galaxy ( M UV )Φ( M UV ) . (36)  \nNote that our galaxy fraction estimates are based on various spectroscopic catalogs, and it is possible that the estimates are biased due to a variety of different spectroscopic selection functions. In order to check this possibility, we compares the estimated galaxy fractions at z ∼ 4 with those derived from the VVDS data (Le F'evre et al. 2013), whose targets are purely selected based on their i -band magnitude. As shown in Figure 9, the estimated galaxy fractions are consistent with each other within the errors, indicating that our estimates are not significantly biased. Future large spectroscopic surveys such as Subaru/Prime Focus Spectrograph (PFS) will allow us to acculately determine the spectroscopic galaxy fractions.  \nFigure 10 and Table 5 show our estimates of the galaxy UV luminosity functions at z ∼ 4 -7. We confirm that  \nour results are consistent with the previous results in the UV magnitude range fainter than -22 mag. This is because the number density of AGNs are negligibly small compared to that of galaxies in this magnitude range. In the brighter magnitude range of M UV < -23 mag, our estimates at z ∼ 4 -7 appear to have bright end excesses of number densities compared to the exponential decline, although the uncertainties are large. The number densities of bright galaxies at z ∼ 6 are determined more precisely than those at z ∼ 5 thanks to the rich spectroscopic data at z ∼ 6 mainly taken in the SHELLQs project. Note that the effect of the Eddington bias (Eddington 1913) should be small on these bright end excesses, because their magnitude ranges are much brighter than the limiting magnitudes, as discussed in Ono et al. (2018).", '4.2.2. Fitting the Galaxy Luminosity Function': "To characterize the derived galaxy UV luminosity functions, we compare our estimates with the following three functions, a DPL function, the Schechter function, and a lensed Schechter function. The forms of the DPL and Schechter functions are already presented in Equations (33) and (35). The lensed Schechter function is a modified Schechter function that considers the effect of gravitational lens magnification by foreground sources (e.g., Wyithe et al. 2011; Takahashi et al. 2011; Mason et al. 2015b; Barone-Nugent et al. 2015). To take into account the magnification effect on the observed shape of the galaxy UV luminosity functions, we basically follow the method presented by Wyithe et al. (2011) and Ono et al. (2018). A gravitationally lensed Schechter function can be estimated from the convolution between the intrinsic Schechter function and the magnification distribution of a Singular Isothermal Sphere (SIS), dP/dµ , weighted by the strong lensing optical depth τ m , which is the fraction of strongly lensed random lines of sight. The overall magnification distribution can be modeled by using the probability distribution function for magnification of multiply lensed sources over a fraction τ m of the sky. To conserve total flux on the cosmic sphere centered on an observer, we need to consider the de-magnification of unlensed sources:  \nµ demag = 1 -〈 µ mult 〉 τ m 1 -τ m , (37)  \nwhere 〈 µ mult 〉 = 4 is the mean magnification of multiply lensed sources. For a given luminosity function φ ( L ), a gravitationally lensed luminosity function φ lensed ( L ) can be obtained by  \nφ lensed ( L ) = (1 -τ m ) 1 µ demag φ ( L µ demag ) + τ m ∫ ∞ 0 dµ 1 µ ( dP m , 1 dµ + dP m , 2 dµ ) φ ( L µ ) , (38)  \nwhere  \ndP m , 1 dµ = { 2 ( µ -1) 3 (for µ > 2) 0 (for 0 < µ < 2) (39)  \nis the magnification distribution as a function of magnification factor µ for the brighter image in a strongly  \nTable 6 Fit Parameters for Luminosity Functions of All Rest-UV Selected Sources (Galaxy+AGN)  \n- ∗ The z ∼ 3 values are based on our fitting for results in Moutard et al. (2020), Bouwens et al. (2021), and Zhang et al. (2021).\n- † The value in parenthesis is fixed to results in Matsuoka et al. (2018c) assuming the redshift evolution of φ ∗ ∝ 10 -0 . 7( z -6) to z = 6 . 8.  \nFigure 8. Evolution of the rest-frame UV luminosity functions of all rest-UV selected sources (galaxies+AGNs) from z ∼ 7 to z ∼ 0. The bottom panel shows the luminosity functions at z ∼ 0 -7, and the black, grey, brown, purple, blue, green, orange, and red symbols show results at z ∼ 0 . 2, 1, 2, 3, 4, 5, 6, and 7, respectively. The circles at z ∼ 4 -7 show our results based on the HSC-SSP survey data, and those at z ∼ 0 -3 are taken from Moutard et al. (2020) at 0 . 05 < z < 0 . 3, 0 . 9 < z < 1 . 3, 1 . 8 < z < 2 . 5, and 2 . 5 < z < 3 . 5. The diamonds are results for AGNs in Zhang et al. (2021), Akiyama et al. (2018), Niida et al. (2020), and Matsuoka et al. (2018c) at z ∼ 3, 4, 5, and 6, respectively, and the squares show results for galaxies in Bouwens et al. (2021) and Ishigaki et al. (2018). The lines show the best-fit DPL+DPL functions at z ∼ 3 -7 (this study), and Schechter functions at z ∼ 0 -2 (Moutard et al. 2020). The top panel shows ratios of the number densities at z ∼ 3 -5 relative to those at z ∼ 6. The shaded regions correspond to the 2 σ uncertainties of the luminosity function parameters. The number density of typical galaxies ( M UV = M ∗ UV ∼ -21 mag) increases only by a factor of ∼ 3 from z ∼ 6 to 3, while the number density of typical quasars ( M UV = M ∗ UV ∼ -27 mag) significantly increases by a factor of ∼ 100.  \n<!-- image -->  \nFigure 9. Spectroscopic galaxy fraction at z ∼ 4. The black circles are those estimated by using the all available spectroscopic catalogs (same as the top panel of Figure 5), and the red circles are those estimated with the VVDS data (Le F'evre et al. 2013).  \n<!-- image -->  \nlensed system given for an SIS and  \ndP m , 2 dµ = 2 ( µ +1) 3 (for µ > 0) (40)  \nis the magnification probability distribution of the second image. Here we consider two cases of results of optical depth estimates to cover a possible range of systematic uncertainties. One is based on the high-resolution raytracing simulations of Takahashi et al. (2011). From their results of the probability distribution function of lensing magnification, the optical depth by foreground sources are estimated to be τ m = 0 . 00231, 0 . 00315, 0 . 00380, and 0 . 00446 at z = 4, 5, 6, and 7, respectively. The other is based on a calibrated Faber-Jackson relation (Faber & Jackson 1976) obtained by Barone-Nugent et al. (2015): τ m = 0 . 0041, 0 . 0054, 0 . 0065, and 0 . 0072 at z = 4, 5, 6, and 7, respectively. Note that these optical depth estimates would be upper limits, because some fraction of lensed dropout sources might be too close to foreground lensing galaxies to be selected as dropouts in our samples. For the Schechter function parameters, we adopt the best-fit values obtained in the Schechter function fitting.  \nIn Figure 10, we show the best-fit functions of these three functional forms with the obtained galaxy UV luminosity function results. Table 7 summarizes the bestfit parameters and the reduced χ 2 values. We find that the DPL and the lensed Schechter functions provide better fits to the observed galaxy UV luminosity functions than the original Schechter functions. The bright-end shapes of the observed galaxy UV luminosity functions cannot be explained by the Schechter functions. The significances of the bright end excess of the number density beyond the Schechter functions are 2 . 9 σ , 1 . 9 σ , 2 . 8 σ , and 2 . 0 σ at z ∼ 4, 5, 6, and 7, respectively. Note that the significances are lower than those in Ono et al. (2018) at z ∼ 4 and 7 because this time we consider the uncertainties of the spectroscopic galaxy fractions, which are not taken into account in Ono et al. (2018). The physical origin of this bright end excess of the number density beyond the Schechter function will be discussed in Section 6.3. The DPL function provides a better fit to the data points than the lensed Schechter function at z ∼ 4 -6, although the significance of this difference is low. The significances of the excess beyond the lensed Schechter functions are 2 . 5(2 . 7) σ , 1 . 4(1 . 6) σ , 2 . 1(2 . 4) σ ,  \nand 1 . 4(1 . 6) σ at z ∼ 4, 5, 6, and 7, respectively, for the optical depth of Barone-Nugent et al. (2015) (Takahashi et al. 2011), slightly smaller than those beyond the Schechter functions.", '4.2.3. Redshift Evolution': 'Figure 11 shows our galaxy luminosity functions at z ∼ 4 -7 with those of Bouwens et al. (2021) at z ∼ 4 -10 and of Bowler et al. (2020) at z ∼ 8 -10. Although Bowler et al. (2020) do not subtract AGN contributions from their estimated luminosity functions, the number densities of their bright sources at z ∼ 8 -10 are likely dominated by galaxies, not by quasars, given the rapid decrease of the quasar luminosity function from z ∼ 3 to 6 as discussed in Section 4.1.4. Figure 11 suggests that the number density of typical galaxies ( M UV = M ∗ UV ∼ -21 mag) significantly increases by a factor of ∼ 100 from ∼ 10 to 4, while that of faint galaxies ( M UV ∼ -16 mag) mildly increases by a factor of ∼ 10. Our comparison also shows that the number density of the bright galaxies at -25 glyph[lessorsimilar] M UV glyph[lessorsimilar] -23 mag does not significantly change at z ∼ 4 -10, and is consistent with no evolution within the 2 σ errors. This is consistent with Bowler et al. (2021), who report little evolution of the number density of bright galaxies at z glyph[greaterorsimilar] 5, although in their comparison they do not subtract AGN contributions from luminosity functions at z ∼ 5 -7. This agreement is expected because Bowler et al. (2021) compare the number densities of M UV > -24 mag sources that are not dominated by AGNs. As shown in Section 4.1.4, the number density of M UV glyph[lessorsimilar] -24 mag sources is dominated by AGNs and evolves rapidly, and we need to subtract AGN contributions to fairly compare the galaxy luminosity functions at the bright end.', '5.1. Angular Correlation Function': 'We calculate angular correlation functions to evaluate the clustering strength of galaxies at z ∼ 2 -6. We use the galaxy samples at z ∼ 1 . 7, 2 . 2, 3, 4, 5, and 6 constructed in Sections 3.1 and 3.5. Note that we do not remove AGNs from our source catalogs. To test the dependence of the clustering strength on the luminosity, we divide our galaxy samples into subsamples by UV magnitude thresholds ( m th UV ). The number of dropouts in the subsamples and their magnitude thresholds are summarized in Table 8. We do not use sources brighter than m cut UV = 20 . 0 mag at each redshift in our analysis. Changing this cut to a fainter magnitude (e.g., 23.0 mag) to remove AGNs does not change results of the angular correlation functions within the errors, because the number of such bright sources is small (see Harikane et al. 2018a). Note that in the calculations we do not use sources in some part of the fields in the Wide layer whose depths are shallow.  \nWe calculate observed angular correlation functions of the subsamples, ω obs ( θ ), using an estimator proposed by Landy & Szalay (1993),  \nω obs ( θ ) = DD ( θ ) -2 DR ( θ ) + RR ( θ ) RR ( θ ) , (41)  \nwhere DD ( θ ), DR ( θ ), and RR ( θ ) are the numbers of galaxy-galaxy, galaxy-random, and random-random  \nFigure 10. Rest-frame UV luminosity functions of galaxies that take into account AGN fraction correction at z ∼ 4, 5, 6, and 7 from top to bottom. The red line is the best-fit DPL function, and the green line is the best-fit Schechter function without considering the lensing effect. The blue shaded region corresponds to the lensed Schechter functions with the two cases of optical depth estimates (Takahashi et al. 2011 and Barone-Nugent et al. 2015). The red circles show our results based on the HSC-SSP survey data. The magenta squares and orange diamonds are results for galaxies taken from Bouwens et al. (2021) and Ishigaki et al. (2018), respectively.  \n<!-- image -->  \nUV [mag]  \nTable 7 Fit Parameters for Galaxy Luminosity Functions  \nFigure 11. Evolution of the rest-frame UV luminosity functions of galaxies from z ∼ 10 to z ∼ 4. The bottom panel shows the luminosity functions at z ∼ 4 -10, and the blue, green, orange, red, magenta, purple, and pink symbols show results at z ∼ 4, 5, 6, 7, 8, 9, and 10, respectively. The circles at z ∼ 4 -7 show our results based on the HSC-SSP survey data, and those at z ∼ 8 -10 are taken from Bowler et al. (2020). The squares show results for galaxies in Bouwens et al. (2021). The lines show the best-fit DPL functions in this study at z ∼ 4 -7 and Bowler et al. (2020) at z ∼ 8 -10. The top panel shows ratios of the number densities relative to those at z ∼ 10. The shaded regions at z ∼ 4 -7 correspond to the 2 σ uncertainties of the luminosity function parameters. The number density of typical galaxies ( M UV = M ∗ UV ∼ -21 mag) significantly increases by a factor of ∼ 100 from ∼ 10 to 4, while that of faint galaxies ( M UV ∼ -16 mag) mildly increases by a factor of ∼ 10. The number density of the bright galaxies at -25 glyph[lessorsimilar] M UV glyph[lessorsimilar] -23 mag does not significantly change at z ∼ 4 -10.  \n<!-- image -->  \npairs normalized by the total numbers of pairs. We use the random catalog whose surface number density is 100 arcmin -2 with the same geometrical shape as the observational data including the mask positions (Coupon et al. 2018). Corrections for contaminations (e.g., Equation (29) in Harikane et al. (2016)) are not applied because the clustering strength of the interlopers is not well-measured. We calculate angular correlation functions in individual fields, and obtain the best-estimate that is the mean weighted by the effective area in each field. Figures 12 and 13 show our calculated angular correlation functions of the subsamples at z ∼ 2 -3 and 4 -6, respectively. We compare our obtained correlation functions with the literature in Figure 14. Our correlations functions are in good agreement with those of Adelberger et al. (2005), Savoy et al. (2011), and Hildebrandt et al. (2009).  \nDue to the finite size of our survey fields, the observed correlation functions is underestimated by a constant value known as the integral constraint, IC (Groth & Peebles 1977). Including a correction for the number of objects in the sample, N (Peebles 1980), the true angular correlation function is given by  \nω ( θ ) = ω obs ( θ ) + IC + 1 N . (42)  \nWe estimate the integral constraint with  \nIC = Σ i RR ( θ i ) ω model ( θ i ) Σ i RR ( θ i ) , (43)  \nwhere ω model ( θ ) is the best-fit model of the correlation function, and i refers the angular bin. IC and ω model ( θ ) are simultaneously determined in the model fitting in Section 5.2.  \nWe estimate statistical errors of the angular correlation functions using the Jackknife estimator. We divide each subsample into Jackknife samples of about 1000 2 arcsec 2 , whose size is larger than the largest angular scale in the correlation function. Removing one Jackknife sample at a time for each realization, we compute the covariance matrix as  \nC ij = N Jack -1 N Jack N Jack ∑ l =1 [ ω l ( θ i ) -¯ ω ( θ i ) ] [ ω l ( θ j ) -¯ ω ( θ j ) ] . (44)  \nwhere N Jack is the total number of the Jackknife samples, ω l is the estimated correlation function from the l th realization, and ¯ ω is the mean correlation function. We apply a correction factor given by Hartlap et al. (2007) to an inverse covariance matrix in order to compensate for the bias introduced by the noise. The inverse of the square root of the inverse covariance matrix is plotted in Figures 12, 13, 14 as uncertainties.', '5.2. Halo Occupation Distribution (HOD) Model Fitting': "We use an HOD model to investigate the relationship between galaxies and their dark matter halos. The HOD model is an analytic framework quantifying a probability distribution of the number of galaxies in dark matter halos (e.g., Seljak 2000; Peacock & Smith 2000; Ma & Fry 2000). The key assumption in the HOD model is that the probability depends only on the halo mass, M h . We can analytically calculate correlation functions and  \nnumber densities from the HOD model. Details of the calculations are presented in Harikane et al. (2016).  \nWe fit our HOD model to the observed angular correlation functions and number densities. In the fitting procedures, the best-fit parameters are determined by minimizing the χ 2 value,  \nχ 2 = ∑ i,j [ ω ( θ i ) -ω model ( θ i )] C -1 i,j [ ω ( θ j ) -ω model ( θ j )] + [ log n obs g -log n model g ] 2 σ 2 log n g , (45)  \nwhere C -1 i,j is the inverse covariance matrix, n g is a space number density of galaxies in the subsample, and σ log n g is its error. We calculate the number density of galaxies corrected for incompleteness using the galaxy UV luminosity functions derived in this work (Section 4.2) and Bouwens et al. (2021). The galaxy number density of each subsample is presented in Table 8. We assume 10% fractional uncertainties in the number densities as Zheng et al. (2007). This 10% uncertainty is a conservative assumption, because the actual statistical uncertainty is typically less than 5%. We constrain the parameters of our HOD model using the Markov Chain Monte Carlo (MCMC) parameter estimation technique.  \nIn our HOD model, an occupation function for central galaxies follows a step function with a smooth transition,  \n〈 N cen ( M h ) 〉 = 1 2 [ 1 + erf ( log M h -log M min √ 2 σ log M h )] . (46)  \nAn occupation function for satellite galaxies is expressed by a power law with a mass cut,  \n〈 N sat ( M h ) 〉 = 〈 N cen ( M h ) 〉 ( M h -M cut M sat ) α . (47)  \nThe total occupation function is  \n〈 N tot ( M h ) 〉 = 〈 N cen ( M h ) 〉 + 〈 N sat ( M h ) 〉 . (48)  \nThese functional forms are motivated by N-body simulations, smoothed particle hydrodynamic simulations, and semi-analytic models for both low and high redshift galaxies (e.g., Kravtsov et al. 2004; Zheng et al. 2005; Garel et al. 2015). Indeed, previous studies demonstrate that this HOD model can explain observed angular correlation functions of high redshift galaxies (Harikane et al. 2016, 2018a; Ishikawa et al. 2017).  \nWe calculate the mean dark matter halo mass of galaxies including both the central and satellite galaxies, 〈 M h 〉 , effective galaxy bias, b eff g , and the satellite fraction, f sat , as follows:  \n〈  \nM h 〉 = 1 n g ∫ dM h dn dM h ( M h , z ) N tot ( M h ) M h , (49) b eff g = 1 n g ∫ dM h dn dM h ( M h , z ) N tot ( M h ) b h ( M h , z ) , (50) f sat = 1 n g ∫ dM h dn dM h ( M h , z ) N sat ( M h ) , (51)  \nwhere dn dM h ( M h , z ), b h ( M h , z ), and n g are the halo mass function, halo bias, and the galaxy number density in the  \nFigure 12. Obtained angular correlation functions of the subsamples and their best-fit by the HOD models at z ∼ 1 . 7, 2 . 2, and 3. The black squares denote the angular correlation function of the subsample at each redshift. The data denoted by the open squares are not used in our HOD model fitting because they are at too-small scales or possibly affected by the non-linear halo bias effect. The dashed and dot-dashed lines represent the one-halo and two-halo terms, respectively, and the solid line is the summations of the one-halo and two-halo terms.  \n<!-- image -->  \nFigure 13. Same as Figure 12 but at z ∼ 4, 5, and 6.  \n<!-- image -->  \nTable 8 Summary of the Clustering Measurements with the HOD Model.  \nNote . -(1) Mean redshift. (2) Threshold apparent magnitude in the rest-frame UV band. (3) Number of galaxies in the subsample. (4) Number density of galaxies in the subsample in units of Mpc -3 . (5) Threshold absolute magnitude in the rest-frame UV band. (6) SFR corresponding to M th UV after the extinction correction in units of M glyph[circledot] yr -1 . (7) Best-fit value of M min in units of M glyph[circledot] . (8) Best-fit value of M sat in units of M glyph[circledot] . The value in parenthesis is derived from M min via Equation (53). (9) Mean halo mass in units of M glyph[circledot] . (10) Effective bias. (11) Satellite fraction. (12) Reduced χ 2 value.  \nmodel (Equation (51) in Harikane et al. 2016), respectively. We use the Behroozi et al. (2013) halo mass function, which is a modification of the Tinker et al. (2008) mass function, and is calibrated at z > 3, the NFW dark matter halo profile (Navarro et al. 1996, 1997), the Duffy et al. (2008) concentration parameter, and the Smith et al. (2003) non-linear matter power spectrum.  \nSome theoretical studies claim that the halo bias is scale-dependent in the quasi-linear scale of r ∼ 50 Mpc (the non-linear halo bias effect; Reed et al. 2009; Jose et al. 2013, 2016, 2017). However, in this study, we assume the scale-independent linear halo bias of Tinker et al. (2010), b ( M h , z ). Instead, we do not use the angular  \ncorrelation functions at 10 '' < θ < 90 '' (10 '' < θ < 120 '' ) in the UltraDeep and Deep (Wide) layers, because they could be affected by the non-linear halo bias effect. We also do not use the measurements at θ ≤ 2 '' that are possibly affected by the source confusion.  \nThe HOD model has five parameters, M min , σ log M h , M cut , M sat , and α . We take M min and M sat as free parameters, which control typical masses of halos having one central and satellite galaxies, respectively, as previous studies (Harikane et al. 2016, 2018a). We fix σ log M h = 0 . 2 and α = 1 . 0, following results of previous studies (e.g., Kravtsov et al. 2004; Zheng et al. 2005; Conroy et al. 2006; Ishikawa et al. 2017). To derive M cut  \nFigure 14. Comparisons of the angular correlation functions with the literature. The black squares represent the correlation functions in this study. The open circles, squares, and diamonds denote results in Adelberger et al. (2005), Savoy et al. (2011), and Hildebrandt et al. (2009), respectively. Our obtained measurements agree well with these previous studies.  \n<!-- image -->  \nFigure 16 shows our results of the halo mass, M h , at z ∼ 2 -6 as a function of the UV magnitude, M UV , with previous studies (Harikane et al. 2016, 2018a) at z ∼ 4 -7. We plot M min and M th UV as M h and M UV , respectively. Table 9 summarizes the results of this work at z ∼ 2 -6 and of Harikane et al. (2016). We find that the new results obtained in this work are consistent with our previous measurements in Harikane et al. (2018a), which are indicated by the crosses in Figure 16. The halo mass of z ∼ 4 -6 galaxies identified in the HSC data ranges from 3 × 10 11 M glyph[circledot] to 1 × 10 13 M glyph[circledot] , which is more massive than those of galaxies identified in the Hubble data (Harikane et al. 2016). The combination of the Hubble and HSC data allows us to investigate the M UV -M h relation over two orders of magnitude in the halo mass at z ∼ 4 and 5. Thus in the following discussion, we will mainly use the results of this work and of Harikane et al. (2016). There is a positive correlation between the UV luminosities and the halo masses at all redshifts, indicating that more UV-luminous galaxies reside in more massive halos, as suggested by previous studies. The slope of the M UV -M h relation becomes steeper at the brighter magnitude, which is similar to the local M ∗ -M h relation (e.g., Leauthaud et al. 2012; Behroozi et al. 2013, 2019; Moster et al. 2013, 2018; Coupon et al. 2015). Note that the uncertainty of the halo mass at z ∼ 6 is as small as those at z ∼ 4 -5 albeit with the large errors in the correlation function measurement at z ∼ 6, because the halo mass is mainly determined by the number density whose uncertainty is 10%, and the slope of the halo mass function is very steep at high redshift.  \n<!-- image -->  \nUV  \nFigure 15. Comparison of the number densities between the HOD models and observations. The gray, brown, purple, blue, green, orange, and red squares (circles) represent the relative differences of the number densities between the HOD models and observations for the subsamples in this work (in Harikane et al. 2016), at z ∼ 1 . 7, 2 . 2, 3, 4, 5, 6, and 7, respectively, as a function of the threshold absolute UV magnitudes, M th UV .  \nfrom M min , we use the relation  \nM cut = M -0 . 5 min , (52)  \nwhich is given by Coupon et al. (2015). Because the exact value of M cut has very little importance compared to the other parameters, this assumption does not change any of our conclusions. For subsamples whose correlation functions are not accurately determined due to the small number of galaxies (subsamples of z ∼ 2 . 2 m th UV = 24 . 0 mag, z ∼ 3 m th UV = 24 . 0 mag, and z ∼ 6 m th UV = 24 . 5 , 25 . 5), we also use the following relation calibrated with results in Harikane et al. (2018a).  \nlog M sat = 3 . 16 log M min -24 . 33 . (53)  \nWe plot the observed angular correlation functions and predictions from their best-fit HOD models in Figures 12 and 13. The best-fit parameters and their 1 σ errors are presented in Table 8. The HOD models can reproduce the observed correlation functions at small (2 glyph[lessorsimilar] θ glyph[lessorsimilar] 10 '' ) and large ( θ glyph[greaterorsimilar] 100 '' ) scales. However, the models underpredict the correlation functions by a factor of 1 . 5 -6 in 10 '' glyph[lessorsimilar] θ glyph[lessorsimilar] 100 '' , the transition scale between 1- and 2-halo terms (the quasi-linear scale), especially in the subsamples at z glyph[greaterorsimilar] 3. These results indicate that the correlation functions at 10 '' glyph[lessorsimilar] θ glyph[lessorsimilar] 100 '' can not be explained by the scale-independent halo bias due to the non-linear halo bias effect in this quasi-linear scale (Jose  \net al. 2013, 2016, 2017). We also find that the best-fit HOD models slightly underpredict the correlation functions of the subsamples of z ∼ 3 and m UV < 25 . 5 and 26 . 0 at θ ≥ 100 '' , although the reduced χ 2 values are not bad ( χ 2 / dof = 1 . 05 , 1 . 54). We have tried to fit the correlation functions of these subsamples with taking M min , σ log M h , M sat , and α as free parameters, but the results does not significantly change. In Figure 15, we compare the observed number densities and predictions from the best-fit HOD models, showing good agreement.", '5.3. M UV -M h Relation': 'We find a redshift evolution of the M UV -M h relation from z ∼ 1 . 7 to 7. For example, M h monotonically decreases from z ∼ 7 to 1 . 7 (from z ∼ 6 to 2 . 2) by a factor of 5 (9) at M UV = -19 . 5 ( -20 . 5). This redshift evolution indicates that the dust-uncorrected SFR increases  \nTable 9 Halo mass and SFR/ ˙ M h in This Study.  \nNote . - (1) Mean redshift. (2) Threshold absolute magnitude in the rest-frame UV band. (3) SFR corresponding to M th UV after the extinction correction in units of M glyph[circledot] yr -1 . (4) Dark matter halo mass ( M min ) in units of M glyph[circledot] . (5) Ratio of the SFR to the dark matter accretion rate in units of 10 -2 .  \nFigure 16. M UV -M h relation. The gray, brown, purple, blue, green, orange, and red filled diamonds (circles) denote the halo masses as a function of the UV magnitude at z ∼ 1 . 7, 2 . 2, 3, 4, 5, 6, and 7, respectively, for the subsamples in this work (in Harikane et al. 2016). We plot M min and M th UV as M h and M UV , respectively. The crosses are results of the previous work based on the early HSC-SSP data (Harikane et al. 2018b). The solid curves show the best-fit relations of Equation (54).  \n<!-- image -->  \nUV  \nwith increasing redshift at fixed dark matter halo mass. We also plot the best-fit M UV -M h relations at z ∼ 1 . 7, 2 . 2, 3, 4, 5, 6, and 7 in Figure 16. These relations are expressed with the following double power law function:  \nM h = M h , 0 2 [ 10 -0 . 4( M UV -M UV , 0 ) α +10 -0 . 4( M UV -M UV , 0 ) β ] (54)  \nwhere M h , 0 and M UV , 0 are characteristic halo mass and UV magnitude, respectively, and α and β are faint and bright end power-law slopes, respectively. In Figure 16, we use parameter sets of (log M h , 0 , M UV , 0 , α, β ) = (11 . 67 , -19 . 30 , 0 . 47 , 2 . 11), (11 . 82 , -19 . 77 , 0 . 61 , 2 . 34), (12 . 05 , -20 . 75 , 0 . 56 , 2 . 15), (11 . 92 , -20 . 90 , 0 . 60 , 2 . 29), (11 . 62 , -20 . 58 , 0 . 51 , 1 . 67), (11 . 44 , -20 . 35 , 0 . 36 , 1 . 54), and (11 . 32 , -20 . 35 , 0 . 36 , 1 . 54) for z ∼ 1 . 7, 2 . 2, 3, 4, 5, 6, and 7, respectively. At z ∼ 7, M UV , 0 , α , and β are fixed to the values at z ∼ 6.  \nIn Figure 17, we compare the mean halo masses, 〈 M h 〉 , of our subsamples with the literature. Because most of the previous studies assume the cosmological parameter set of (Ω m , Ω Λ , h, σ 8 ) = (0 . 3 , 0 . 7 , 0 . 7 , 0 . 9) that is different from our assumption, we obtain HOD model fitting results for our data with (Ω m , Ω Λ , h, σ 8 ) = (0 . 3 , 0 . 7 , 0 . 7 , 0 . 9) for comparison. Similarly, the results of the previous studies are re-calculated with the same cosmological parameter sets if different cosmological parameter set is assumed. In this way, we conduct our comparisons using an equivalent set of cosmological parameters across all datasets. In Figure 17, we find that our results at z ∼ 3 and 4 are consistent with those of the previous studies within the uncertainties. While the  \nFigure 17. Comparison of the mean dark matter halo masses, 〈 M h 〉 , with the literature under the same cosmology. We plot M th UV as M UV . The diamonds represent the mean dark matter halo masses in this work with the cosmological parameters of ( h, Ω m , Ω Λ , σ 8 ) = (0 . 7 , 0 . 3 , 0 . 7 , 0 . 9). The circles are the results of Harikane et al. (2016) with the same cosmology. The black symbols denote results of the previous studies. We plot the results of Hildebrandt et al. (2009, squares), Lee et al. (2006, upward triangles), Ouchi et al. (2005, stars), and Hamana et al. (2004, downward triangle). The downward triangles have no error bars, because Hamana et al. (2004) do not provide errors of the mean dark matter halo mass. We also show the results of Barone-Nugent et al. (2014) as black open squares, which are re-calcutaed with the cosmological parameters of ( h, Ω m , Ω L , σ 8 ) = (0 . 7 , 0 . 3 , 0 . 7 , 0 . 9). The mean dark matter halo mass of the faintest subsample at z ∼ 3 is slightly more massive than that of the next faintest subsample, because of the higher fraction of satellite galaxies that are typically reside in massive halos (see Table 8).  \n<!-- image -->  \nprevious results at z ∼ 5 are largely scattered, our z ∼ 5 results are placed near the center of the distribution of the previous studies. At z ∼ 6, our results agree with that of Barone-Nugent et al. (2014). In summary, our results are consistent with most of the previous studies. Furthermore, our results improve on both the statistics and the dynamic range covered in M UV .', '5.4. SFR/ ˙ M h -M h Relation': 'We estimate a ratio of the SFR to the dark matter accretion rate, SFR/ ˙ M h , or the baryon conversion efficiency, SFR/ ( f b ˙ M h ), where f b = Ω b / Ω m is the cosmic baryon fraction. Since the baryon gas accretes into the halo together with dark matter, this ratio indicates the star formation efficiency. In this paper, the star formation efficiency indicates SFR/ ˙ M h or SFR/ ( f b ˙ M h ), not the ratio of the SFR to the gas mass ( SFR/M gas ), which is usually used in radio astronomy. We derive the dustuncorrected SFRs ( SFR UV ) from UV luminosities using the calibration used in Madau & Dickinson (2014) with the Salpeter (1955) IMF:  \nSFR UV ( M glyph[circledot] yr -1 ) = 1 . 15 × 10 -28 L UV (erg s -1 Hz -1 ) . (55)  \nWe correct the SFR for the dust extinction using an attenuation-UV slope ( β UV ) relation (Meurer et al. 1999) and β UV -M UV relation at each redshift. We use the β UV -M UV relation in Bouwens et al. (2014) at z glyph[greaterorsimilar] 4 and linearly extrapolate the relation with fixing the slope at z glyph[lessorsimilar] 4. The estimated SFRs are presented in Table 8. We calculate ˙ M h as a function of halo mass and redshift using an analytic formula obtained from N-body simulation results in Behroozi & Silk (2015, Equation (B8)). Note that the accretion rates in Behroozi & Silk (2015) are typically ∼ 2 times lower than those calculated based on Equations (E2)-(E6) in Behroozi et al. (2013), which are used in our previous work (Harikane et al. 2018a), because the Behroozi et al. (2013) accretion rates only trace the progenitors of z = 0 halos.  \nWe plot the dust-corrected SFR/ ˙ M h ratios at z ∼ 2 -7 as a function of the halo mass in Figure 18. The results  \nFigure 18. SFR / ˙ M h and baryon conversion efficiency (SFR / ( f b ˙ M h )) as a function of the halo mass. The gray, brown, purple, blue, green, orange, and red filled diamonds (circles) denote the ratios as a function of the halo mass, at z ∼ 1 . 7, 2 . 2, 3, 4, 5, 6, and 7, respectively, for the subsamples in this work (in Harikane et al. 2016). The statistical errors for our data are smaller than the symbols (diamonds). The black solid curve is the fitting formulae of Equation (56) for the SFR/ ˙ M h -M h relation at z ∼ 2 -7, and the gray shaded region represents the 2 σ typical scatter (0 . 3 dex) of the data points compared to the relation.  \n<!-- image -->  \n/circledot  \nare also summarized in Table 9. The black solid curve in Figure 18 represents the following SFR/ ˙ M h -M h relation:  \nSFR ˙ M h = 2 × 3 . 2 × 10 -2 ( M h / 10 11 . 5 ) -1 . 2 +( M h / 10 11 . 5 ) 0 . 5 (56)  \nThis relation agrees with the measured SFR/ ˙ M h ratios at z ∼ 2 -7 within 0.3 dex that is a typical 2 σ scatter. This good agreement indicates that the star formation efficiency does not significantly change beyond 0.3 dex in the wide redshift range of z ∼ 2 -7, suggesting the existence of the fundamental relation between the star formation and the mass accretion (the growth of the galaxy and its dark matter halo assembly), as discussed in Harikane et al. (2018a) at z ∼ 4 -7 (see also Bian et al. 2013).  \nFigure 19. Comparison of the SFR / ˙ M h ratio with the literature. The red solid curve shows the redshift-dependent relation constrained in this work (Equation (57)), and the shaded region represents the 2 σ typical uncertainties. Results compared include those from our previous work (Equation (25) in Harikane et al. 2018a, red dashed curve), Behroozi et al. (2019, blue solid curve), Moster et al. (2018, gray dashed curve), Tacchella et al. (2018, green solid curve), Mason et al. (2015a, gray dotted curve), and Behroozi et al. (2013, blue dashed curve). All results are converted to use the Salpeter (1955) IMF (equation (55)). The result of Harikane et al. (2018a) is re-calculated based on the accretion rate of Behroozi & Silk (2015), and thus is ∼ 2 times higher than the original relation in Harikane et al. (2018a).  \n<!-- image -->  \n/circledot  \n/circledot  \n/circledot  \nOn the other hand, the SFR/ ˙ M h ratio gradually increases with increasing redshift within 0.3 dex from z ∼ 5 to 1 . 7. If we take this possible evolution into account, the ratio can be expressed as,  \nSFR ˙ M h = 2 × 3 . 2 × 10 -2 ( M h / 10 11 . 5 ) -1 . 2 +( M h / 10 11 . 5 ) 0 . 5 × (0 . 53 tanh[0 . 54(2 . 9 -z )] + 1 . 53) (57)  \nThis relation indicates that the star formation efficiency does not significantly change from z ∼ 7 to 5, and then gradually increases within a factor of ∼ 2 from z ∼ 5 to 1 . 7, still consistent with the results of Harikane et al. (2018a), who have identified redshift-independent relation at z ∼ 4 -7 within 0.15 dex. The reason for this elevated efficiency at z < 5 is not clear. One possibility is an increase of the metallicity in galaxies toward lower redshift, resulting in more efficient gas cooling.  \nIn Figure 19, we compare our SFR / ˙ M h -M h relation with results in the literature (Behroozi et al. 2013; Mason et al. 2015a; Harikane et al. 2018a; Tacchella et al. 2018; Moster et al. 2018; Behroozi et al. 2019). The relations in the literature show similar trends to our result; the SFR / ˙ M h ratio has a peak of SFR / ˙ M h ∼ 0 . 1 -0 . 01 around the halo mass of 10 11 -10 12 M glyph[circledot] . However, the slopes of the relation at the high-mass and low-mass ends are different between these studies. These differences are possibly due to differences in used observational datasets, the halo mass functions, and details of modeling. We also note that there are several systematic uncertainties on the SFR / ˙ M h measurements. Assumptions on the IMF and dust attenuation have impacts on the SFR. The conversion factor between the SFR and UV luminosity (Equation (55)) depends on the stellar age and metallicity (e.g., Madau & Dickinson 2014).', '6.1. Physical Origin of the Cosmic SFR Density Evolution': 'In Section 5.4, we find the fundamental SFR/ ˙ M h -M h relation; the value of SFR/ ˙ M h at fixed M h does not significantly change beyond 0 . 3 dex at z ∼ 2 -7. We examine whether this fundamental SFR/ ˙ M h -M h relation is consistent with the observational results, i.e.,  \ncosmic SFR densities and the UV luminosity functions. We calculate the cosmic SFR density as follows:  \nρ SFR = ∫ dM h dn dM h SFR = ∫ dM h dn dM h ˙ M h SFR ˙ M h , (58)  \nwhere SFR/ ˙ M h at z ∼ 2 -7 is obtained as a function of M h in Section 5.4 (Equations (56) or (57)). We integrate down to the halo mass corresponding to the SFR of 0 . 3 M glyph[circledot] yr -1 ( M UV = -17 mag with the Madau & Dickinson 2014 calibration), as previous studies (Bouwens et al. 2015, 2020; Finkelstein et al. 2015b; Oesch et al. 2018).  \nFirstly we assume the redshift-independent SFR/ ˙ M h -M h relation (Equation (56)). Figure 20 compares our calculated SFR densities (the dashed curve) with observational results in the literature (Madau & Dickinson 2014; Finkelstein et al. 2015b; McLeod et al. 2016; Oesch et al. 2018; Bouwens et al. 2020). The results in the literature are all converted to use the calibration of Madau & Dickinson (2014) with the Salpeter (1955) IMF (Equation (55)). We find that our calculation well reproduces the overall trend of the cosmic SFR density evolution; the calculated density increases from z ∼ 10 to 4 -2, and decreases from z ∼ 4 -2 to 0. However, the SFR densities are underpredicted compared to the observations at z ∼ 1 -2 by ∼ 0 . 3 dex.  \nThen we use the gradually evolving SFR/ ˙ M h -M h relation (Equation (57)), instead of Equation (56). As shown in Figure 20, our calculated cosmic SFR densities (the solid curve) based on Equation (57) agree well with the observations especially at z ∼ 1 -2, compared to the calculation based on the redshift-independent relation (the dashed curve). Quantitatively, the reduced χ 2 value improves significantly from χ 2 / dof = 25 . 4 to 3 . 0. These analyses indicate that the overall trend of the redshift evolution can be reproduced by the redshift-independent SFR/ ˙ M h -M h relation (Equation (56)), but the gradual increase of the star formation efficiency at z < 5 (Equation (57)) is needed to quantitatively reproduce the observed SFR densities. Note that we have constrained the SFR/ ˙ M h -M h relations by using normal star forming galaxies (dropout, BX, and BM galaxies), and not considered quiescent or dusty starburst galaxies. Quiescent  \nFigure 21. Mechanism of the cosmic SFR density evolution. Top panel: the purple, blue, green, orange, and red curves indicate the number density of halos whose masses are M h = 10 10 , 10 11 , 10 12 , 10 13 , and 10 14 M glyph[circledot] , respectively. The number density is calculated by using the Behroozi et al. (2013) halo mass function. The black curve represents a weighted number density based on Equation (57). Middle panel: same as the top panel but for the dark matter accretion rate calculated by using the formula in Behroozi & Silk (2015). Bottom panel: same as the top panel but for the cosmic SFR density. The gray circles are observed cosmic SFR densities taken from Madau & Dickinson (2014), Finkelstein et al. (2015b), McLeod et al. (2016), and Bouwens et al. (2020). Since the cosmic SFR density at a given halo mass is proportional to the galaxy (halo) number densities and SFRs (accretion rates) as shown in Equation (58), the redshift evolution of the cosmic SFR density is made by the monotonic steep increase of the halo number density from z ∼ 10 to z ∼ 4 and the monotonic decrease of the accretion rate from z ∼ 2 to z ∼ 0, resulting a peak of the SFR density around z ∼ 2 -3.  \n<!-- image -->  \nz  \nFigure 20. Cosmic SFR density. The bottom panel shows the comparison of the cosmic SFR densities. The red curves with the shade represent cosmic SFR densities derived from Equation (58). For the solid and dashed lines, we assume the redshift-dependent and independent SFR/ ˙ M h -M h relations (Equations (57) and (56)), respectively, which are constrained in this paper at 2 glyph[lessorsimilar] z glyph[lessorsimilar] 7. The shaded regions correspond to the 1 σ (0.15 dex) scatter in the SFR/ ˙ M h -M h relation. We integrate down to the halo mass corresponding to the SFR of 0 . 3 M glyph[circledot] yr -1 ( M UV = -17 mag with the Madau & Dickinson 2014 calibration), as previous studies (Bouwens et al. 2015, 2020; Finkelstein et al. 2015b; Oesch et al. 2018). The black circles show observed cosmic SFR densities taken from Madau & Dickinson (2014), Finkelstein et al. (2015b), McLeod et al. (2016), and Bouwens et al. (2020). All results are converted to use the calibration of Madau & Dickinson (2014) with the Salpeter (1955) IMF (Equation (55)). The top panel shows ∆ SFR ˙ M h = SFR ˙ M h ( z ) / SFR ˙ M h ( z = 5) at a given halo mass. The sold and dashed curves are calculated with Equations (57) and (56), respectively.  \nand dusty starburst galaxies are expected to have lower and higher SFR/ ˙ M h ratios than the normal star forming galaxies, and effects of these galaxies are thought to be less dominant at z glyph[greaterorsimilar] 4 (e.g., Muzzin et al. 2013; Bouwens et al. 2020; Fudamoto et al. 2021). Although quiescent and dusty starburst galaxies would be non-negligible at z glyph[lessorsimilar] 3, the good agreement at z ∼ 1 -2 still indicates that the star formation efficiency averaged over all galaxy populations gradually increases at z < 5, as long as the observed SFR densities (including the treatment of the dust extinction correction) are correct.  \nThe good agreement with the overall trend indicates that the evolution of the cosmic SFR densities is primarily driven by the monotonic increase of the halo number density and the monotonic decrease of the accretion rate, given the weak redshift evolution of the SFR/ ˙ M h -M h relation, as discussed in Harikane et al. (2018a). The number density of halos at a given halo mass increases due to structure formation from z ∼ 10 to a certain redshift at z glyph[lessorsimilar] 4 depending on the mass and then becomes almost constant after that (the top panel in Figure 21), resulting in the increase of the galaxy number density from z ∼ 10 to z glyph[lessorsimilar] 4. The dark matter (and gas) accretion rate monotonically decreases over the whole redshift  \n<!-- image -->  \nz  \nrange due to the cosmic expansion, with a steep drop from z ∼ 2 to 0 (the middle panel in Figure 21), resulting in the monotonic decrease SFR of each galaxy at a given halo mass. Because the cosmic SFR density at a given halo mass is proportional to the number density and mass accretion rate (or SFR) as shown in Equation (58), the calculated cosmic SFR density has a peak at z ∼ 2 -3 (the bottom panel in Figure 21). More specifically, the product of the number density and mass accretion rate for each halo mass has a peak at a cer-  \n0  \ntain redshift due to the increase of the number density and decrease of the accretion rate, with the peak redshift depending on the halo mass, and the SFR/ ˙ M glyph[circledot] -M h relation determines the peak redshift of the cosmic SFR density integrated over the halo mass, as shown in the bottom panel in Figure 21).  \nIn our calculation we integrate Equation (58) down to the SFR of 0 . 3 M glyph[circledot] yr -1 , corresponding to the halo mass of ∼ 3 × 10 10 M glyph[circledot] at z ∼ 7. This integration limit is chosen to match the calculations at z glyph[greaterorsimilar] 2 in previous studies (Bouwens et al. 2015, 2020; Finkelstein et al. 2015b; Oesch et al. 2018), and slightly different from the calculations in Madau & Dickinson (2014), who integrate down to 0 . 03 L ∗ . This difference does not affect our discussions above. Bouwens et al. (2021) show that the M ∗ UV parameter of the luminosity function is almost constant ( M ∗ UV ∼ -21 mag) from z ∼ 8 to 2, and then decreases to M ∗ UV ∼ -18 mag toward z ∼ 0. At z ∼ 2 -8, the SFR corresponding to 0 . 03 L ∗ with M ∗ UV = -21 mag is roughly 0 . 5 M glyph[circledot] yr -1 , comparable to our integration limit. Below z ∼ 2, the corresponding SFR is smaller than our limit, e.g., 0 . 03 M glyph[circledot] yr -1 with M ∗ UV ∼ -18 mag at z ∼ 0. However, even if we integrate Equation (58) down to this SFR limit, the calculated cosmic SFR density increases only by 0.1 dex at z < 3, and cannot explain the 0.3 dex difference between the observed SFR densities and the dashed curve in Figure 20.  \nWe also calculate the UV luminosity function at each redshift as follows:  \nΦ( M UV ) = dn dM h dM h dM UV . (59)  \nFrom Equation (57), we can obtain the M UV -M h relation and dM h dM UV at each redshift, since ˙ M h can be expressed as a function of M h and z (Behroozi & Silk 2015). Note that M UV is the observed absolute magnitude after dust extinction assuming the attenuation-UV slope ( β UV ) relation (Meurer et al. 1999) and β UV -M UV relations (Bouwens et al. 2014). We correct for satellite galaxies using satellite fractions measured in previous studies (Wake et al. 2011; Martinez-Manso et al. 2015; McCracken et al. 2015; Harikane et al. 2016, 2018a), although the correction is not large at high redshift. The calculated UV luminosity functions at z ∼ 0 -10 are plotted in Figure 22. We find that the calculated luminosity functions are in rough agreement with observed results given the 0 . 15 dex (1 σ ) uncertainty in SFR/ ˙ M h , indicating that our SFR/ ˙ M h -M h relation is consistent with the observed redshift evolution of the UV luminosity function.  \nMason et al. (2015a) and Tacchella et al. (2018) also report that the constant star formation efficiency model can reproduce the UV luminosity functions at z glyph[greaterorsimilar] 4. This is consistent with our results of the SFR/ ˙ M h -M h relation in this study and our previous work (Harikane et al. 2018a). Bouwens et al. (2021) claim that the evolution of the luminosity function at z ∼ 2 . 5 -10 can be explained by the halo mass function and the constant star formation efficiency model. This is qualitatively consistent with our results, but the gradual increase of the star formation efficiency at z < 5 is needed to quantitatively reproduce the redshift evolution of the cosmic SFR  \nFigure 22. Comparison of the rest-frame UV luminosity functions. Solid curves are the calculated luminosity functions from Equation (59) with the SFR / ˙ M h -M h relation of Equation (57). The points show the observed galaxy luminosity functions of Bouwens et al. (2021, circles), Moutard et al. (2020, squares), and Arnouts et al. (2005, stars). The shaded regions correspond to the 1 σ (0.15 dex) scatter in the SFR/ ˙ M h -M h relation.  \n<!-- image -->  \ndensity, as discussed above.', '6.2. Future Prospects for Star Formation at z > 10': 'In the previous section, we show that the constant star formation efficiency (Equation (56)) can reproduce the evolution of the cosmic SFR density at 5 glyph[lessorsimilar] z glyph[lessorsimilar] 10. By assuming that this SFR/ ˙ M h -M h relation does not evolve to higher redshifts, we can predict the cosmic SFR density at z > 10 based on the evolution of the halo mass function and dark matter accretion rate.  \nFigure 23 compares our calculated SFR densities (the red curve) with predictions from models in the literature (Mason et al. 2015a, Mashian et al. 2016, Sun & Furlanetto 2016, Tacchella et al. 2018, Oesch et al. 2018, and Behroozi et al. 2020). We find that the cosmic SFR density based on the constant star formation efficiency rapidly decreases with increasing redshift as ∝ 10 -0 . 5(1+ z ) , similar to the predictions of other models. More quantitatively, the SFR densities from observations at z glyph[lessorsimilar] 10 and our predictions at z glyph[greaterorsimilar] 10 are well-fitted with the following function:  \nρ SFR / [M glyph[circledot] yr -1 Mpc -3 ] = 1 61 . 7 × (1 + z ) -3 . 13 +1 . 0 × 10 0 . 22(1+ z ) +2 . 4 × 10 0 . 50(1+ z ) -3 . 0 , (60)  \nas shown in Figure 24. This is contrast to the extrapolation of the fitting function in Madau & Dickinson (2014) that shows a smooth decline as ∝ (1 + z ) -2 . 9 at z > 10 (the gray dashed curve in Figures 23 and 24), and possible estimates of the SFR densities at z > 12 based on z ∼ 6 passive galaxies in Mawatari et al. (2020). James Webb Space Telescope (JWST) will allow us to directly observe galaxies at z ∼ 10 -20, and investigate whether the SFR density rapidly decreases as predicted in many models or not, providing insights into star formation efficiencies in z ∼ 10 -20 galaxies.  \nUsing this method to predict the SFR density at z > 10, we can obtain a rough estimate of the epoch of the  \n<!-- image -->  \nz  \nFigure 23. Comparison of the cosmic SFR density at z > 7. The red curve with the shade represents the cosmic SFR density calculated in this work based on the constant star formation efficiency at z > 5 (Equation (57)), integrated down to the SFR of 0 . 3 M glyph[circledot] yr -1 ( M UV = -17 mag), as previous studies (Bouwens et al. 2015, 2020; Finkelstein et al. 2015b; Oesch et al. 2018). The gray dashed curve shows the extrapolation of the relation of Madau & Dickinson (2014) at z > 6. The other curves show predictions from models of Mason et al. (2015a, blue dashed curve), Mashian et al. (2016, orange dashed curve), Sun & Furlanetto (2016, green solid curve), Tacchella et al. (2018, green dashed curve), and Behroozi et al. (2020, blue solid curve). The orange shaded region indicates a prediction of the halo evolution model in Oesch et al. (2018). All results are converted to use the Salpeter (1955) IMF (equation (55)).Figure 24. Fit to the observed cosmic SFR densities. The red dashed curve represents our fit (Equation (60)) to the observed cosmic SFR densities at z glyph[lessorsimilar] 10 and the calculated SFR densities at z > 10 in this work. The gray dashed curve shows the fit in Madau & Dickinson (2014). All results are converted to use the Salpeter (1955) IMF (equation (55)).  \n<!-- image -->  \nfirst star formation. We calculate the cumulative number of formed stars as a function of redshift:  \nN star ( z ) = V survey ∫ t ( z ) t ( z ) -t LF dtρ SFR /M FS , (61)  \nwhere V survey is the survey volume, t LF is the typical lifetime of the star, and M FS is a typical mass of the first star. Here we assume t LF = 3 Myr (Schaerer 2002) and M FS = 100 M glyph[circledot] (e.g., Hirano et al. 2015). We adopt V survey = (3 h -1 Mpc) 3 that is the volume of the simulation box in Hirano et al. (2015). We calculate ρ SFR using Equation (58) by extrapolating the SFR/ ˙ M h -M h re-  \ntion both to the higher redshift and lower mass range. We integrate down to the halo mass of 10 5 M glyph[circledot] , comparable to halo masses of first stars in simulations (e.g., Hirano et al. 2015). Figure 25 shows the calculated cumulative number of formed stars. The shaded region indicates possible uncertainties of the low-mass slope of the SFR/ ˙ M h -M h relation and the mass limit of the integration that we adopt 10 4 -10 6 M glyph[circledot] . The number reaches 1 around z ∼ 16 -27, implying the first star formation in this epoch. This formation epoch agrees with theoretical simulations (Hirano et al. 2014, 2015), although this is a very rough estimate. Especially, it is not clear whether the assumed SFR/ ˙ M h -M h relation holds at z > 10 or not, because physics in the first star formation are expected to be different from star/galaxy formation at the epoch we currently observe due to the evolution of physical parameters such as metallicity.', '6.3. Origin of the Bright End Excess of the Galaxy Luminosity Function': "As presented in Section 4.2.2, the obtained galaxy UV luminosity functions cannot be explained by the Schechter functions at the bright end ( M UV glyph[lessorsimilar] -23 mag), indicating the existence of the bright end excess of the number density beyond the Schechter function. Since these luminosity functions are estimated based on the spectroscopic galaxy fractions, the bright end excess cannot be explained by apparent AGNs (e.g., quasars). Here we discuss the following five possibilities for the origin of the bright end excess, (A) lensed galaxies, (B) mergers, (C) inefficient mass quenching, (D) low dust obscuration and (E) hidden AGN activity.  \n(A) Lensed galaxies. In Section 4.2.2, the lensed Schechter function can provide better fits than the  \nFigure 25. Cumulative number of formed stars as a function of the redshift. The red line shows the cumulative number of formed stars calculated based on the SFR density (Equation (61)). We extrapolate the SFR/ ˙ M h -M h relation both to the higher redshift and lower mass range, and the shaded region indicates possible uncertainties of the SFR/ ˙ M h ratio and the mass limit of the integration. The cumulative number reaches 1 around z ∼ 16 -27, implying the first star formation at this epoch.  \n<!-- image -->  \nSchechter functions. This indicates that the bright end excess can be explained by gravitationally lensed galaxies that are apparently bright. However, in our visual inspection of galaxies located at the bright end, we do not find a significant number of lensed galaxies that show elongated morphologies. In addition, the z ∼ 4 luminosity function still shows the bright end excess at the 2 . 5 -2 . 7 σ levels beyond the lensed Schechter function, similar to the results of Ono et al. (2018). As discussed in Section 4.2.2, the calculated lensed Schechter function is an upper limit, because some fraction of lensed galaxies might be too close to foreground lensing galaxies to be selected as dropouts in our samples. Thus the bright end excess does not seem to be easily explained by lensed galaxies, although high spatial resolution images are required to draw definitive conclusions.  \n- (B) Mergers. Some studies suggest that the major merger at high redshift is an important phase for formation of massive galaxies (e.g., Sawicki et al. 2020). Due to the limited spatial resolutions of the HSC images, some merging galaxies are not resolved and look like one bright galaxy. If most of galaxies at the bright end turn out to be merging galaxies, the bright end excess can be explained by merging galaxies. High spatial resolution images such as those obtained by Hubble are indeed useful to investigate this possibility (Bowler et al. 2017), but only a small fraction of our sources are observed with Hubble . Instead, Shibuya et al. (2021) have made high spatial resolution images by using the super-resolution technique, and estimate a major merger fraction of bright galaxies in galaxy samples of Ono et al. (2018). They find that the major merger fraction is 10 -70% for bright galaxies ( -24 glyph[lessorsimilar] M UV glyph[lessorsimilar] -22 mag) at z ∼ 4 -7. This  \ncan partly explain the number density excess of relatively bright galaxies ( M UV glyph[greaterorsimilar] -23 . 6 mag), but is not sufficient to explain the excess of very bright galaxies ( M UV glyph[lessorsimilar] -23 . 6 mag). Thus, it seems that the bright end excess cannot be explained by mergers.  \n- (C) Inefficient mass quenching. It is thought that the exponential cutoff at the bright end of the Schechter function is caused by the mass quenching that suppresses star formation activity in massive halos (e.g., AGN feedback, virial shock heating). Thus the bright end excess beyond the Schechter function implies inefficient mass quenching in these high redshift bright galaxies as discussed in Ono et al. (2018) and Bowler et al. (2020). To investigate this possibility, we obtain the star formation efficiency, SFR/ ˙ M h , as a function of the halo mass in cases of the Schechter and DPL luminosity functions. Using the halo mass function and galaxy luminosity function, we obtain the M UV -M h relation by the abundance matching technique:  \n∫ ∞ M h dM ' h dn dM h ( M ' h , z )(1 + f sat ) = ∫ M UV -∞ dM ' UV Φ galaxy ( M ' UV , z ) . (62)  \nThen we calculate the SFR and mass accretion rate in the same manner as Section 5.4, and obtain the SFR/ ˙ M h ratio.  \nFigure 26 presents the luminosity function, the M UV -M h relation, and the star formation efficiency (the SFR/ ˙ M h -M h relation) in cases of the DPL and Schechter functions at z ∼ 4, 6, 8, and 10. We find that estimated halo masses from our abundance matching agree well with those from the clustering analysis (the middle panel in Figure 26), as discussed in Harikane et al. (2016). If we assume the steep decline of the star formation efficiency ( SFR/ ˙ M h ) toward the massive end (the dashed line in the right panel), the calculated number density shows the exponential decline at the bright end similar to the Schechter function (the dashed line in the left panel), and cannot reproduce the excess of the observed number densities. To reproduce the bright end excess of the number density like the DPL function (the solid line in the left panel), higher star formation efficiencies are needed at the massive end (the solid line in the right panel), compared to the case of the Schechter function. For example, ∼ 2 times higher star formation efficiency is needed in halos of M h glyph[similarequal] 10 13 (10 12 ) M glyph[circledot] at z ∼ 6 (8). The high star formation efficiency at the bright end can be made by the inefficient mass quenching. Indeed the mass quenching is expected to be less efficient at higher redshift because of the shorter time scale of gas cooling and/or weaker AGN feedback due to the decreasing number of AGNs (i.e., quasar luminosity function) as discussed in Section 4.1.4.  \n- (D) Low dust obscuration. Bowler et al. (2020) dis-  \nFigure 26. Galaxy luminosity function (left), M UV -M h relation (middle), and star formation efficiency ( SFR/ ˙ M h -M h relation, right) in cases of the DPL and Schechter luminosity functions. The blue, orange, magenta, and pink solid (dashed) lines represent relations in the case of the DPL (Schechter) luminosity functions at z ∼ 4, 6, 8, and 10, respectively. To reproduce the observed bright end excess of the number density like the DPL function (the solid line in the left panel), the higher star formation efficiency (the solid line in the right panel) is needed at the massive end, compared to the case of the Schechter function (the dashed line). The blue and orange symbols in the middle panel show the results of this study (diamonds), Harikane et al. (2018a, crosses), and Harikane et al. (2016, circles). In the right panel, the black solid curve is the fitting formulae of Equation (56), and the gray shaded region represents the 2 σ typical scatter (0 . 3 dex).  \n<!-- image -->  \ncuss the possibility that the intrinsic (without dust attenuation) UV luminosity function has a shallower decline at the bright end, and the dust obscuration controls the shape of the luminosity function. In the calculations above (C), we assume the attenuationβ UV relation in Meurer et al. (1999) and the β UV -M UV relations in Bouwens et al. (2014). However, the attenuation curve of high redshift galaxies is not well-understood (e.g., Hashimoto et al. 2019; Fudamoto et al. 2020; Harikane et al. 2020b; Bakx et al. 2020), and the β UV -M UV relation is not well-constrained at this very bright magnitude range (i.e., ∼ -23 mag). In addition, some recent studies report dust-poor UV luminous star forming galaxies at z > 2 (MarquesChaves et al. 2020, 2021). Thus it is possible that the dust obscuration in these bright galaxies is lower than what we assumed, resulting the bright end excess beyond the Schechter function.  \n- (E) Hidden AGN activity. Although we subtract the number density of AGNs by the spectroscopic galaxy fractions, it is possible that there are still hidden AGNs in the galaxy luminosity function. UV luminosities of such hidden AGNs could be boosted due to AGN activity (e.g., Kim & Im 2021), resulting in the bright end excess. AGN activity of such sources can be probed only by deep spectroscopy covering several high ionization lines. Indeed, some studies report possible AGN activity in bright ( M UV glyph[lessorsimilar] -22 mag) galaxies at z glyph[greaterorsimilar] 7 (e.g., Laporte et al. 2017; Mainali et al. 2018; Jiang et al. 2021; Endsley et al. 2021; Onoue et al. 2021). Since we do not know the fraction of such hidden AGNs in our sample due to the lack of deep spectroscopic data at the bright end, we cannot rule out the possibility that hidden AGNs make the bright end excess.  \nBased on these discussions above, we conclude that the bright end excess is possibly made by (C) inefficient mass quenching, (D) low dust obscuration and/or (E) hidden AGN activity, although it is possible that the dominant effect in making the bright end excess changes with red-  \n/circledot  \nshift. We cannot distinguish these possibilities with the current datasets, and future large and deep observations are needed. For example, Euclid and Nancy Grace Roman Space Telescope can identify a large number of high redshift galaxies located at the bright end, allowing us to investigate the clustering and the star formation efficiency of such bright galaxies. ALMA follow-up observations for a statistical sample of bright galaxies will reveal the typical dust properties of these galaxies. Deep spectroscopy for a large number of bright galaxies with Subaru/PFS will allow us to investigate hidden AGN activity in such bright galaxies.", '7. Summary': "In this paper, we have identified 1,978,462 dropout candidates at z ∼ 4 -7 from ∼ 300 deg 2 deep optical imaging data obtained in the HSC-SSP survey. Among these dropout candidates, 1037 dropouts are spectroscopically identified in our follow-up observations and the literature. Typical contamination rates are < 20% and < 40% at z ∼ 4 -5 and z ∼ 6 -7, respectively. Combined with z ∼ 2 -3 galaxy samples, we have a total of 4,100,221 sources at z ∼ 2 -7, which is the largest sample of high redshift galaxies to date. Using this sample, we have calculated the luminosity functions and angular correlation functions, and investigated statistical properties of these sources.  \nOur major findings are summarized below:  \n- 1. We have obtained rest-frame UV luminosity functions at z ∼ 4 -7 (Figures 5 and 6). Combined with results based on the complementary ultra-deep Hubble data and wide-area SDSS data, we have probed the luminosity function in a very wide UV magnitude range of -29 glyph[lessorsimilar] M UV glyph[lessorsimilar] -14 mag, corresponding to the luminosity range of 0 . 002 L ∗ UV glyph[lessorsimilar] L UV glyph[lessorsimilar] 2000 L ∗ UV .\n- 2. Spectroscopic galaxy fractions indicate that most of the sources are AGNs (galaxies) at M UV < -24 ( M UV > -22) mag (Figures 5 and 6). We have found that the luminosity function in this very wide magnitude range can be well-fitted by the DPL+DPL or DPL+Schechter functions (Figure  \n- 7), indicating that the dropout luminosity function is a superposition of the AGN luminosity function (dominant at the bright end) and the galaxy luminosity function (dominant at the faint end).\n- 3. We have estimated the galaxy luminosity functions by subtracting the AGN contributions using the spectroscopic galaxy fractions. The obtained galaxy luminosity functions show the bright end excess of the number density beyond the Schechter function at glyph[greaterorsimilar] 2 σ levels (Figure 10), which is possibly made by the inefficient mass quenching, low dust obscuration, and/or hidden AGN activity (Section 6.3).\n- 4. We have derived angular correlation functions of galaxies at z ∼ 2 -6 (Figures 12 and 13). Combined with the HOD model analyses and previous clustering measurements for faint galaxies at z ∼ 4 -7, we have obtained the relation between the dark matter halo mass and the UV magnitude over two orders of magnitude in the halo mass (Figure 16).\n- 5. We have calculated the ratio of the SFR to the dark matter accretion rate, SFR/ ˙ M h , and identified an SFR/ ˙ M h -M h relation which does not show strong redshift evolution beyond 0.3 dex at z ∼ 2 -7 (Figure 18). This weak evolution indicates that the star formation efficiency does not significantly change at high redshift, and star formation activities are regulated by the dark matter mass assembly, as suggested by our earlier work at z ∼ 4 -7 (Harikane et al. 2018a). Meanwhile, the SFR/ ˙ M h ratio gradually increases with decreasing redshift from z ∼ 5 to 2 within 0.3 dex.\n- 6. We have found that the SFR/ ˙ M h -M h relation can reproduce the redshift evolution of the cosmic SFR density and the UV luminosity function (Figures 20 and 22). These good agreements indicate that the evolution of the cosmic SFR densities is primarily driven by the steep increase of the halo number density from z ∼ 10 to z glyph[lessorsimilar] 4 due to the structure formation, and the decrease of the accretion rate due to the cosmic expansion with a steep drop from z ∼ 2 to 0 (Figure 21).  \nWe have further showed that the cosmic SFR density at z > 10 decreases towards higher redshift more rapidly than the extrapolation of the fitting function in Madau &Dickinson (2014) if we assume the constant star formation efficiency (Figures 23 and 24). JWST observations allow us to directly investigate this rapid evolution by measuring the cosmic SFR densities at z > 10, providing insights into star formation efficiency in the early universe.  \nWe thank the anonymous referee for a careful reading and valuable comments that improved the clarity of the paper. We thank Masao Hayashi, Yoshiyuki Inoue, Kentaro Nagamine, Rui Marques Coelho Chaves, Shun Saito, Masayuki Tanaka, Yoshihiko Yamada for useful comments and discussions, and Thibaud Moutard and Yechi  \nZhang for providing their data of UV luminosity functions. This work is supported by the World Premier International Research Center Initiative (WPI Initiative), the Ministry of Education, Culture, Sports, Science and Technology, Japan, as well as KAKENHI Grant-in-Aid for Scientific Research (19J01222, 20H00180, 21H04467, and 21K13953) through the Japan Society for the Promotion of Science (JSPS).  \nThe HSC collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from the Japanese Cabinet Office, MEXT, the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.  \nThis paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at http://dm.lsst.org  \nThis paper is based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center (ADC) at NAOJ. Data analysis was in part carried out with the cooperation of Center for Computational Astrophysics (CfCA), NAOJ. This work was supported by the joint research program of the Institute for Cosmic Ray Research (ICRR), University of Tokyo.  \nPart of the data were obtained and processed as part of the CFHT Large Area U-band Deep Survey (CLAUDS), which is a collaboration between astronomers from Canada, France, and China described in Sawicki et al. (2019). CLAUDS is based on observations obtained with MegaPrime/ MegaCam, a joint project of CFHT and CEA/DAPNIA, at the CFHT which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. CLAUDS uses data obtained in part through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories, Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance of China. CLAUDS uses data products from TERAPIX and the Canadian Astronomy Data Centre (CADC) and was carried out using resources from Compute Canada and Canadian Advanced Network For Astrophysical Research (CANFAR). CADC is operated by the National Research Council of Canada with the support of the Canadian Space Agency.  \nThe Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics,  \nGarching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lo-  \nrand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.  \nThe Cosmic Dawn Center is funded by the Danish National Research Foundation under grant No. 140. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 847523 'INTERACTIONS'.  \nSoftware: galaxev (Bruzual & Charlot 2003), GalSim (v1.4; Rowe et al. 2015), hscpipe (Bosch et al. 2018), MIZUKI (Tanaka 2015; Tanaka et al. 2018), SynPipe (Huang et al. 2018)  \nAPPENDIX  \nCatalog of Spectroscopic Sources  \nTable 10 summarizes 1037 spectroscopically confirmed galaxies and AGNs in our dropout samples at z ∼ 4 -7.  \nTable 10  \nSpectroscopically Identified Galaxies and AGNs in Our Dropout Samples.  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nTable 10 -Continued  \nNote . - (1) Object ID. (2) Right ascension. (3) Declination. (4) Spectroscopic redshift. (5) Apparent magnitude. (6) Absolute UV magnitude. (7) The dropout sample in which the source is selected: 1 = g -dropout, 2 = r -dropout, 3 = i -dropout, and 4 = z -dropout. (8) Galaxy/AGN flag (1 = galaxy; 2 = AGN). (9) Reference for the spectroscopic redshift.", 'REFERENCES': "Abazajian, K., Adelman-McCarthy, J. K., Agueros, M. A., et al.  \n2004, AJ, 128, 502 Adams, N. J., Bowler, R. A. A., Jarvis, M. J., et al. 2020, MNRAS, 494, 1771 Adelberger, K. L., & Steidel, C. C. 2000, ApJ, 544, 218 Adelberger, K. L., Steidel, C. C., Pettini, M., et al. 2005, ApJ, 619, 697 Adelberger, K. L., Steidel, C. C., Shapley, A. E., et al. 2004, ApJ, 607, 226 Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4 Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, PASJ, 71, 114 Akiyama, M., He, W., Ikeda, H., et al. 2018, PASJ, 70, S34 Alavi, A., Siana, B., Richard, J., et al. 2016, ApJ, 832, 56 Arnouts, S., Schiminovich, D., Ilbert, O., et al. 2005, ApJ, 619, L43 Atek, H., Richard, J., Kneib, J.-P., & Schaerer, D. 2018, MNRAS, 479, 5184 Atek, H., Richard, J., Jauzac, M., et al. 2015, ApJ, 814, 69 Axelrod, T., Kantor, J., Lupton, R. H., & Pierfederici, F. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7740, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 15  \n- Ba˜nados, E., Venemans, B. P., Decarli, R., et al. 2016, ApJS, 227, 11  \nBakx, T. J. L. C., Tamura, Y., Hashimoto, T., et al. 2020, MNRAS, 493, 4294  \n- Barone-Nugent, R. L., Wyithe, J. S. B., Trenti, M., et al. 2015, MNRAS, 450, 1224\n- Barone-Nugent, R. L., Trenti, M., Wyithe, J. S. B., et al. 2014, ApJ, 793, 17\n- Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143\n- Behroozi, P., Conroy, C., Wechsler, R. H., et al. 2020, MNRAS, 499, 5702\n- Behroozi, P. S., & Silk, J. 2015, ApJ, 799, 32\n- Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57 Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38\n- Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393\n- Bian, F., Fan, X., Jiang, L., et al. 2013, ApJ, 774, 28\n- Binney, J. 1977, ApJ, 215, 483\n- -. 2004, MNRAS, 347, 1093\n- Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5 Bouwens, R., Gonz'alez-L'opez, J., Aravena, M., et al. 2020, ApJ, 902, 112\n- Bouwens, R. J., Oesch, P. A., Illingworth, G. D., Ellis, R. S., & Stefanon, M. 2017, ApJ, 843, 129\n- Bouwens, R. J., Stefanon, M., Oesch, P. A., et al. 2019, ApJ, 880, 25  \n- Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2014, ApJ, 793, 115  \n-. 2015, ApJ, 803, 34  \n- Bouwens, R. J., Oesch, P. A., Stefanon, M., et al. 2021, arXiv e-prints, arXiv:2102.07775\n- Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645\n- Bowler, R. A. A., Adams, N. J., Jarvis, M. J., & Haußler, B. 2021, MNRAS, 502, 662\n- Bowler, R. A. A., Dunlop, J. S., McLure, R. J., & McLeod, D. J. 2017, MNRAS, 466, 3612\n- Bowler, R. A. A., Jarvis, M. J., Dunlop, J. S., et al. 2020, MNRAS, 493, 2059\n- Bowler, R. A. A., Dunlop, J. S., McLure, R. J., et al. 2014, MNRAS, 440, 2810  \n-. 2015, MNRAS, 452, 1817  \n- Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000  \nCalzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682  \nCapak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99  \nCastellano, M., Yue, B., Ferrara, A., et al. 2016, ApJ, 823, L40  \n- Cheema, G. K., Sawicki, M., Arcila-Osejo, L., et al. 2020, MNRAS, 494, 804\n- Childress, M. J., Lidman, C., Davis, T. M., et al. 2017, MNRAS, 472, 273\n- Coleman, G. D., Wu, C. C., & Weedman, D. W. 1980, ApJS, 43, 393\n- Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201\n- Cooper, M. C., Griffith, R. L., Newman, J. A., et al. 2012, MNRAS, 419, 3018\n- Coupon, J., Czakon, N., Bosch, J., et al. 2018, PASJ, 70, S7 Coupon, J., Arnouts, S., van Waerbeke, L., et al. 2015, MNRAS, 449, 1352\n- Cowley, W. I., Caputi, K. I., Deshmukh, S., et al. 2018, ApJ, 853, 69\n- Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11\n- Cuby, J.-G., Le F'evre, O., McCracken, H., et al. 2003, A&A, 405, L19\n- Curtis-Lake, E., McLure, R. J., Pearce, H. J., et al. 2012, MNRAS, 422, 1425\n- Dressler, A., Bigelow, B., Hare, T., et al. 2011, PASP, 123, 288 Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008,  \nMNRAS, 390, L64  \n- Eddington, A. S. 1913, MNRAS, 73, 359\n- Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7 Endsley, R., Stark, D. P., Charlot, S., et al. 2021, MNRAS, 502, 6044  \nFaber, S. M., & Jackson, R. E. 1976, ApJ, 204, 668 Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, 1657-1669  \n- Finkelstein, S. L., Song, M., Behroozi, P., et al. 2015a, ApJ, 814, 95\n- Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015b, ApJ, 810, 71\n- Finkelstein, S. L., Bagley, M., Song, M., et al. 2021, arXiv e-prints, arXiv:2106.13813\n- Fudamoto, Y., Oesch, P. A., Faisst, A., et al. 2020, A&A, 643, A4 Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, arXiv e-prints, arXiv:2109.10378  \nFurusawa, H., Kosugi, G., Akiyama, M., et al. 2008, ApJS, 176, 1  \nFurusawa, H., Koike, M., Takata, T., et al. 2018, PASJ, 70, S3  \n- Garel, T., Blaizot, J., Guiderdoni, B., et al. 2015, MNRAS, 450, 1279\n- Garilli, B., McLure, R., Pentericci, L., et al. 2021, A&A, 647, A150  \nGehrels, N. 1986, ApJ, 303, 336  \n- Giallongo, E., Grazian, A., Fiore, F., et al. 2015, A&A, 578, A83 Giavalisco, M. 2002, ARA&A, 40, 579\n- Glikman, E., Djorgovski, S. G., Stern, D., et al. 2011, ApJ, 728, L26\n- Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580\n- Groth, E. J., & Peebles, P. J. E. 1977, ApJ, 217, 385\n- Gunn, J. E., & Stryker, L. L. 1983, ApJS, 52, 121\n- Hamana, T., Ouchi, M., Shimasaku, K., Kayo, I., & Suto, Y. 2004, MNRAS, 347, 813\n- Harikane, Y., Laporte, N., Ellis, R. S., & Matsuoka, Y. 2020a, ApJ, 902, 117\n- Harikane, Y., Ouchi, M., Ono, Y., et al. 2016, ApJ, 821, 123 -. 2018a, PASJ, 70, S11\n- Harikane, Y., Ouchi, M., Shibuya, T., et al. 2018b, ApJ, 859, 84\n- Harikane, Y., Ouchi, M., Ono, Y., et al. 2019, ApJ, 883, 142\n- Harikane, Y., Ouchi, M., Inoue, A. K., et al. 2020b, ApJ, 896, 93  \nHartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399  \n- Hashimoto, T., Inoue, A. K., Mawatari, K., et al. 2019, PASJ, 71, 71\n- Hasinger, G., Capak, P., Salvato, M., et al. 2018, ApJ, 858, 77 Hatfield, P. W., Bowler, R. A. A., Jarvis, M. J., & Hale, C. L. 2018, MNRAS, 477, 3760\n- Hayashi, M., Shimakawa, R., Tanaka, M., et al. 2020, PASJ, 72, 86  \nHiguchi, R., Ouchi, M., Ono, Y., et al. 2019, ApJ, 879, 28  \nHildebrandt, H., Pielorz, J., Erben, T., et al. 2009, A&A, 498, 725  \nHirano, S., Hosokawa, T., Yoshida, N., Omukai, K., & Yorke,  \nH. W. 2015, MNRAS, 448, 568  \nHirano, S., Hosokawa, T., Yoshida, N., et al. 2014, ApJ, 781, 60 Hu, E. M., Cowie, L. L., Songaila, A., et al. 2016, ApJ, 825, L7 Hu, W., Wang, J., Zheng, Z.-Y., et al. 2017, ApJ, 845, L16 Huang, S., Leauthaud, A., Murata, R., et al. 2018, PASJ, 70, S6 Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS,  \n- 442, 1805\n- Inoue, A. K., Hasegawa, K., Ishiyama, T., et al. 2018, PASJ, 70, 55\n- Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12 -. 2018, ApJ, 854, 73\n- Ishikawa, S., Kashikawa, N., Toshikawa, J., et al. 2017, ApJ, 841, 8\n- Ishikawa, S., Kashikawa, N., Tanaka, M., et al. 2020, ApJ, 904, 128\n- Ivezic, Z., Tyson, J. A., Abel, B., et al. 2008, ArXiv e-prints, arXiv:0805.2366\n- Jiang, L., McGreer, I. D., Fan, X., et al. 2016, ApJ, 833, 222\n- Jiang, L., Shen, Y., Bian, F., et al. 2017, ApJ, 846, 134 Jiang, L., Kashikawa, N., Wang, S., et al. 2021, Nature\n- Astronomy, 5, 256\n- Jose, C., Baugh, C. M., Lacey, C. G., & Subramanian, K. 2017, MNRAS, 469, 4428\n- Jose, C., Lacey, C. G., & Baugh, C. M. 2016, MNRAS, 463, 270 Jose, C., Subramanian, K., Srianand, R., & Samui, S. 2013, MNRAS, 429, 2333\n- Kakuma, R., Ouchi, M., Harikane, Y., et al. 2019, arXiv e-prints, arXiv:1906.00173\n- Kashikawa, N., Aoki, K., Asai, R., et al. 2002, PASJ, 54, 819 Kashikawa, N., Ishizaki, Y., Willott, C. J., et al. 2015, ApJ, 798, 28\n- Kawanomoto, S., Uraguchi, F., Komiyama, Y., et al. 2018, PASJ, 70, 66\n- Kim, Y., & Im, M. 2021, ApJ, 910, L11\n- Kim, Y., Im, M., Jeon, Y., et al. 2020, ApJ, 904, 111  \nKnapp, G. R., Leggett, S. K., Fan, X., et al. 2004, AJ, 127, 3553  \nKomiyama, Y., Obuchi, Y., Nakaya, H., et al. 2018, PASJ, 70, S2  \nKonno, A., Ouchi, M., Shibuya, T., et al. 2018, PASJ, 70, S16  \nKravtsov, A. V., Berlind, A. A., Wechsler, R. H., et al. 2004,  \nApJ, 609, 35  \n- Kriek, M., Shapley, A. E., Reddy, N. A., et al. 2015, ApJS, 218, 15\n- Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24 Landy, S. D., & Szalay, A. S. 1993, ApJ, 412, 64\n- Laporte, N., Nakajima, K., Ellis, R. S., et al. 2017, ApJ, 851, 40 Le F'evre, O., Cassata, P., Cucciati, O., et al. 2013, A&A, 559, A14\n- Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219\n- Leauthaud, A., Tinker, J., Bundy, K., et al. 2012, ApJ, 744, 159 Lee, K.-S., Giavalisco, M., Gnedin, O. Y., et al. 2006, ApJ, 642, 63 Lewis, I. J., Cannon, R. D., Taylor, K., et al. 2002, MNRAS, 333, 279\n- Lidman, C., Tucker, B. E., Davis, T. M., et al. 2020, MNRAS, 496, 19  \n- Lilly, S. J., Le Fevre, O., Hammer, F., & Crampton, D. 1996, ApJ, 460, L1\n- Ma, C.-P., & Fry, J. N. 2000, ApJ, 543, 503\n- Madau, P. 1995, ApJ, 441, 18\n- Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415\n- Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388\n- Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20\n- Mainali, R., Zitrin, A., Stark, D. P., et al. 2018, MNRAS, 479, 1180\n- Mallery, R. P., Mobasher, B., Capak, P., et al. 2012, ApJ, 760, 128\n- Marques-Chaves, R., Schaerer, D., ' Alvarez-M'arquez, J., et al. 2021, MNRAS, 507, 524\n- Marques-Chaves, R., ' Alvarez-M'arquez, J., Colina, L., et al. 2020, MNRAS, 499, L105\n- Martin, D. C., Seibert, M., Buat, V., et al. 2005, ApJ, 619, L59 Martinez-Manso, J., Gonzalez, A. H., Ashby, M. L. N., et al. 2015, MNRAS, 446, 169\n- Mashian, N., Oesch, P. A., & Loeb, A. 2016, MNRAS, 455, 2101 Mason, C. A., Trenti, M., & Treu, T. 2015a, ApJ, 813, 21 Mason, C. A., Treu, T., Schmidt, K. B., et al. 2015b, ApJ, 805, 79 Masters, D., Capak, P., Salvato, M., et al. 2012, ApJ, 755, 169 Masters, D. C., Stern, D. K., Cohen, J. G., et al. 2017, ApJ, 841, 111\n- Matsuoka, Y., Onoue, M., Kashikawa, N., et al. 2016, ApJ, 828, 26\n- -. 2018a, PASJ, 70, S35\n- Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2018b, ApJS, 237, 5 Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018c, ApJ, 869, 150\n- Matsuoka, Y., Iwasawa, K., Onoue, M., et al. 2019, ApJ, 883, 183 Mawatari, K., Inoue, A. K., Hashimoto, T., et al. 2020, ApJ, 889, 137\n- McCracken, H. J., Wolk, M., Colombi, S., et al. 2015, MNRAS, 449, 901\n- McLeod, D. J., McLure, R. J., & Dunlop, J. S. 2016, MNRAS, 459, 3812\n- McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696\n- Mehta, V., Scarlata, C., Rafelski, M., et al. 2017, ApJ, 838, 29 Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64 Miyatake, H., Harikane, Y., Ouchi, M., et al. 2021, arXiv e-prints, arXiv:2103.15862\n- Miyazaki, S., Komiyama, Y., Nakaya, H., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 0\n- Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1\n- Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27\n- More, S., Miyatake, H., Mandelbaum, R., et al. 2015, ApJ, 806, 2 Morishita, T., Trenti, M., Stiavelli, M., et al. 2018, ApJ, 867, 150 Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121\n- Moster, B. P., Naab, T., & White, S. D. M. 2018, MNRAS, 477, 1822\n- Moutard, T., Sawicki, M., Arnouts, S., et al. 2020, MNRAS, 494, 1894\n- Murayama, T., Taniguchi, Y., Scoville, N. Z., et al. 2007, ApJS, 172, 523\n- Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18\n- Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563  \n-. 1997, ApJ, 490, 493  \n- Newman, J. A., Cooper, M. C., Davis, M., et al. 2013, ApJS, 208, 5\n- Niida, M., Nagao, T., Ikeda, H., et al. 2020, ApJ, 904, 89 Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104\n- Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labb'e, I., & Stefanon, M. 2018, ApJ, 855, 105\n- Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010, ApJ, 725, L150\n- Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713  \nOno, Y., Ouchi, M., Harikane, Y., et al. 2018, PASJ, 70, S10  \n- Ono, Y., Itoh, R., Shibuya, T., et al. 2021, ApJ, 911, 78\n- Onoue, M., Kashikawa, N., Willott, C. J., et al. 2017, ApJ, 847, L15\n- Onoue, M., Matsuoka, Y., Kashikawa, N., et al. 2021, arXiv e-prints, arXiv:2106.13807\n- Ouchi, M., Shimasaku, K., Okamura, S., et al. 2004, ApJ, 611, 685 Ouchi, M., Hamana, T., Shimasaku, K., et al. 2005, ApJ, 635, L117\n- Ouchi, M., Shimasaku, K., Akiyama, M., et al. 2008, ApJS, 176, 301\n- Ouchi, M., Mobasher, B., Shimasaku, K., et al. 2009, ApJ, 706, 1136\n- Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, ApJ, 723, 869\n- Ouchi, M., Harikane, Y., Shibuya, T., et al. 2018, PASJ, 70, S13 Pˆaris, I., Petitjean, P., Aubourg, ' E., et al. 2018, A&A, 613, A51 Parsa, S., Dunlop, J. S., & McLure, R. J. 2018, MNRAS, 474, 2904\n- Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194\n- Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144 Peebles, P. J. E. 1980, The large-scale structure of the universe (Princeton University Press)\n- Pell'o, R., Hudelot, P., Laporte, N., et al. 2018, A&A, 620, A51 Pentericci, L., Vanzella, E., Castellano, M., et al. 2018, A&A, 619, A147\n- Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13\n- Press, W. H., & Schechter, P. 1974, ApJ, 187, 425\n- Qiu, Y., Wyithe, J. S. B., Oesch, P. A., et al. 2018, MNRAS, 481, 4885\n- Reed, D. S., Bower, R., Frenk, C. S., Jenkins, A., & Theuns, T. 2009, MNRAS, 394, 624\n- Rees, M. J., & Ostriker, J. P. 1977, MNRAS, 179, 541\n- Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015,  \nAstronomy and Computing, 10, 121  \n- Saito, T., Shimasaku, K., Okamura, S., et al. 2008, ApJ, 675, 1076 Salpeter, E. E. 1955, ApJ, 121, 161\n- Savoy, J., Sawicki, M., Thompson, D., & Sato, T. 2011, ApJ, 737, 92\n- Sawicki, M., Arcila-Osejo, L., Golob, A., et al. 2020, MNRAS, 494, 1366\n- Sawicki, M., Arnouts, S., Huang, J., et al. 2019, MNRAS, 489, 5202\n- Sawicki, M. J., Lin, H., & Yee, H. K. C. 1997, AJ, 113, 1\n- Scannapieco, E., & Oh, S. P. 2004, ApJ, 608, 62\n- Schaerer, D. 2002, A&A, 382, 28\n- Schechter, P. 1976, ApJ, 203, 297\n- Schenker, M. A., Robertson, B. E., Ellis, R. S., et al. 2013, ApJ, 768, 196\n- Schlafly, E. F., Finkbeiner, D. P., Juri'c, M., et al. 2012, ApJ, 756, 158\n- Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525\n- Seljak, U. 2000, MNRAS, 318, 203\n- Sharp, R., Saunders, W., Smith, G., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6269, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. I. S. McLean & M. Iye, 62690G\n- Shibuya, T., Miura, N., Iwadate, K., et al. 2021, arXiv e-prints, arXiv:2106.03728  \nShibuya, T., Ouchi, M., & Harikane, Y. 2015, ApJS, 219, 15 Shibuya, T., Ouchi, M., Konno, A., et al. 2018a, PASJ, 70, S14 Shibuya, T., Ouchi, M., Harikane, Y., et al. 2018b, PASJ, 70, S15 Shioya, Y., Taniguchi, Y., Sasaki, S. S., et al. 2009, ApJ, 696, 546 Silk, J. 1977, ApJ, 211, 638  \n- Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311\n- Stefanon, M., Labb'e, I., Bouwens, R. J., et al. 2017, ApJ, 851, 43 -. 2019, ApJ, 883, 99  \nSteidel, C. C., Adelberger, K. L., Giavalisco, M., Dickinson, M., & Pettini, M. 1999, ApJ, 519, 1  \n- Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2003, ApJ, 592, 728\n- Steidel, C. C., Giavalisco, M., Pettini, M., Dickinson, M., & Adelberger, K. L. 1996, ApJ, 462, L17\n- Steidel, C. C., Shapley, A. E., Pettini, M., et al. 2004, ApJ, 604, 534\n- Stevans, M. L., Finkelstein, S. L., Wold, I., et al. 2018, ApJ, 863, 63\n- Sun, G., & Furlanetto, S. R. 2016, MNRAS, 460, 417\n- Tacchella, S., Bose, S., Conroy, C., Eisenstein, D. J., & Johnson, B. D. 2018, ApJ, 868, 92\n- Takahashi, R., Oguri, M., Sato, M., & Hamana, T. 2011, ApJ, 742, 15\n- Tanaka, M. 2015, ApJ, 801, 20\n- Tanaka, M., Coupon, J., Hsieh, B.-C., et al. 2018, PASJ, 70, S9 Tasca, L. A. M., Le F'evre, O., Ribeiro, B., et al. 2017, A&A, 600, A110\n- Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 Tinker, J. L., Robertson, B. E., Kravtsov, A. V., et al. 2010, ApJ, 724, 878\n- Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99  \n- Toshikawa, J., Kashikawa, N., Overzier, R., et al. 2016, ApJ, 826, 114\n- Toshikawa, J., Uchiyama, H., Kashikawa, N., et al. 2018, PASJ, 70, S12\n- van der Burg, R. F. J., Hildebrandt, H., & Erben, T. 2010, A&A,\n- 523, A74 Wake, D. A., Whitaker, K. E., Labb'e, I., et al. 2011, ApJ, 728, 46 Wang, B., & Heckman, T. M. 1996, ApJ, 457, 645 Wang, F., Wu, X.-B., Fan, X., et al. 2016, ApJ, 819, 24 Willott, C. J., Delorme, P., Reyl'e, C., et al. 2009, AJ, 137, 3541 -. 2010, AJ, 139, 906 Willott, C. J., McLure, R. J., Hibon, P., et al. 2013, AJ, 145, 4 Wyithe, J. S. B., Yan, H., Windhorst, R. A., & Mao, S. 2011, Nature, 469, 181 Yang, J., Fan, X., Wu, X.-B., et al. 2017, AJ, 153, 184 Yuan, F., Lidman, C., Davis, T. M., et al. 2015, MNRAS, 452, 3047 Zhang, H., Ouchi, M., Itoh, R., et al. 2020, ApJ, 891, 177 Zhang, Y., Ouchi, M., Gebhardt, K., et al. 2021, arXiv e-prints, arXiv:2105.11497 Zheng, Z., Coil, A. L., & Zehavi, I. 2007, ApJ, 667, 760 Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791"}