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

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-2000L*uv combined with previous studies, revealing that the dropout luminosity function is a superposition of the AGN luminosity function dominant at Muv<-24 mag and the galaxy luminosity function dominant at Muv>-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>2sigma levels, which is 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/(dMh/dt), 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 redshift evolution of the cosmic SFR density, suggesting that the 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 $\propto10^{-0.5(1+z)}$, which will be directly tested with JWST.


Introduction
Studying statistical properties of galaxies is important to understanding 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 correlation 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 formation and feedback processes.
Great 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 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;McLure et al. 2013;Schenker et al. 2013;Bouwens et al. 2015Bouwens et al. , 2019Bouwens et al. , 2021Finkelstein 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., Atek et al. 2015Atek et al. , 2018Ishigaki et al. 2015Ishigaki et al. , 2018Oesch et al. 2015Oesch et al. , 2018Alavi et al. 2016;Castellano 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).
Investigating 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;Bower et al. 2006Croton et al. 2006, inefficiency of gas cooling in massive dark matter halos due to virial shock heating (e.g., Binney 1977;Rees & Ostriker 1977;Silk 1977;Benson et al. 2003), and dust obscuration that 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. 2014Bowler et al. , 2015Bowler et al. , 2017Bowler et al. , 2020Stefanon et al. 2017Stefanon et al. , 2019Morishita et al. 2018;Ono et al. 2018;Stevans 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).
Together with studying luminosity functions, the clustering analysis with the angular correlation function is important to understanding 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;Coupon et al. 2015;More et al. 2015), the abundance matching/empirical model (e.g., Behroozi et al. 2013Behroozi et al. , 2019Moster et al. 2013Moster et al. , 2018Finkelstein et al. 2015a), and the clustering analysis (e.g., Harikane et al. 2016Harikane et al. , 2018aIshikawa et al. 2017Ishikawa et al. , 2020Cowley et al. 2018;Qiu et al. 2018;Cheema et al. 2020). Since the weak-lensing analysis cannot be applied at z  2 due to the limited number of the background galaxies and their lower-image quality with the current observational data sets (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. 2004Ouchi et al. , 2005Hildebrandt et al. 2009;Savoy et al. 2011;Bian et al. 2013; Barone-Nugent et al. 2014;Harikane et al. 2016Harikane et al. , 2018bHatfield et al. 2018). These studies reveal that the more UV-luminous galaxies reside in more massive halos.
Recently, Harikane et al. (2018b) have identified a tight relation between the ratio of the SFR to the dark matter accretion rate,  M SFR 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 redshiftindependent 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  M M SFR h h relation can reproduce the UV luminosity functions at z  4 (e.g., Mason et al. 2015a;Tacchella et al. 2018;Harikane et al. 2018a;Bouwens et al. 2021) and the trend of the redshift evolution of the cosmic SFR density (e.g., Mason et al. 2015a;Oesch et al. 2018;Tacchella et al. 2018;Harikane et al. 2018a), 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  M M SFR h 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.
In this work, we present new measurements of the rest-frame 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 (HSC) Subaru Strategic Program (SSP) survey (Aihara et al. 2018; see also Miyazaki et al. 2012Miyazaki et al. , 2018Furusawa et al. 2018;Komiyama 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 HSC-SSP survey data. One program is our luminous Lyman break galaxy (LBG) or dropout galaxy studies, named Great Optically Luminous Dropout Research Using Subaru HSC (GOLDRUSH; Ono et al. 2018;Toshikawa et al. 2018;Harikane et al. 2018a). The other program is highredshift Lyα emitter (LAE) studies using HSC narrowband filters, named Systematic Identification of LAEs for Visible Exploration and Reionization Research Using Subaru HSC (SILVERRUSH; Inoue et al. 2018;Konno et al. 2018;Ouchi et al. 2018;Harikane et al. 2018bHarikane et al. , 2019Shibuya et al. 2018aShibuya et al. , 2018bHiguchi et al. 2019;Kakuma et al. 2021;Ono et al. 2021). Our new LBG catalogs are made public on our project webpage 23 or Zenodo. 24 This paper is organized as follows. We show the observational data sets 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).

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. 2019). 25 The HSC-SSP survey obtains deep optical imaging data with the five broadband filters, g, r, i, z, and y , 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 ). The applied flags and threshold values are summarized in Table 2. In addition to these flags, we mask some regions that are affected by the bright source halos or the ghosts of bright sources not flagged.
The HSC data are reduced by the HSC-SSP collaboration with hscPipe ) that is the HSC data reduction pipeline based on the Large Synoptic Survey Telescope pipeline (Axelrod et al. 2010;Ivezic et al. 2019). 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 (Schlafly et al. 2012;Tonry et al. 2012;Magnier et al. 2013). PSFs are calculated in hscPipe, and typical full width at half maximum of the PSFs are 0″.6-0″.9.
We 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. 2018;Toshikawa et al. 2018;Harikane et al. 2018a), 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 that 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 that 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. We 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 idropouts, 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 images 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.

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. 2003Steidel et al. , 2004Adelberger 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 with 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.

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 (Lewis et al. 2002;Sharp et al. 2006) on the Anglo-Australian Telescope from 2018 December 31 to 2019 January 3 (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). False grizy None of the pixels in the footprint of an object is labeled as bad mask_brightstar_any False grizy None of the pixels in the footprint of an object is close to bright sources mask_brightstar_ghost15 −99 grizy None of the pixels in the footprint of an object is close to the ghost masks sdsscentroid_flag False ri for g-drop Object centroid measurement has no problem False True ri for g-drop Detected in r and i. False/True g/iz for r-drop Undetected in g and detected in r and i. False/True gr/zy for i-drop Undetected in g and r, and detected in z and y. False/True gri/y for z-drop Undetected in g, r, and i, and detected in y. blendedness_abs_flux < 0.2 ri for g-drop The target photometry is not significantly affected by neighbors.
The number of exposures is equal to or larger than 3/5 in gr/izy.
In the DEIMOS observations, we used the 600ZD grating with the GG455 filter. The spectroscopic observations 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 s 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 were 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 Å. Flux calibration was achieved with data of a standard star G191B2B. The details of the DEIMOS observations will be presented in Y. Ono et al. 2021, in preparation. In the AAOmega+2dF observations, we used the X5700 dichroic beam splitter, the 580V grating with the central wavelength at 4821 Å in the blue channel, and the 385R grating with the central wavelength at 7251 Å in the red channel. This configuration covered a wavelength range of 3800-8800 Å 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 s 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).
In 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 s per each mask. The FOCAS data were reduced with the focasred pipeline. Wavelength calibrations were 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 Å. Flux calibration was performed with data of a standard star BD+28d4211. The details of the FOCAS observations will be presented in Y. Ono et al. 2021, in preparation. A 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.
In addition to the observations described above, we include results of our observations with the Inamori Magellan Areal  11-14, 2008November 29-30, December 1-2, December 18-20, 2009October 11-13, 2010 February 8-9, July 9-10, and 2011 January 3-4. In these observations, the main targets were high-redshift LAE candidates found in the deep Subaru Suprime-Cam narrowband images obtained in the SXDS (Ouchi et al. 2008(Ouchi et al. , 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 Spectroscopy (COSMOS) pipeline. 26

Source Selection at z ∼ 4-7
From 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.
First, to identify secure sources, we select sources whose signal-to-noise ratios (S/Ns) are higher than 5 within 1″.5 diameter apertures in the i band for g-dropouts and in the z 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 e = 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 dotted-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). 26 http://code.obs.carnegiescience.edu/cosmos band for r-dropouts. For i-dropouts, we select sources with S/Ns higher than 5 and 4 in the z and y bands, respectively, because the y-band images are relatively shallow. For the zdropouts, we select sources with S/Ns 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;Toshikawa et al. 2018;Harikane et al. 2018a), we adopt the following color criteria:g- -> z y 1.6. 10 ( ) As shown in Figure 1, these color criteria are set to avoid lowredshift 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 the previous studies similarly using the z−y color (e.g., Ouchi et al. 2009).
To remove foreground interlopers, we exclude sources with continuum detections at >2σ levels in the g band for rdropouts, 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.
Using 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, 2,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  4) galaxy population to date. In particular, 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 the 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.

Spectroscopically Identified Sources
In our samples, a total of 46 sources are identified as z > 3 objects through our spectroscopic follow-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   Å, confirmed in our FOCAS observations. This source is 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α In addition, we incorporate results of our spectroscopic observations for high-redshift galaxies with Magellan/IMACS. We also check spectroscopic catalogs in other studies Figure 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 Å bins with red filled circles to show the Lyman break features. (Cuby et al. 2003;Ouchi et al. 2008;Saito et al. 2008;Willott et al. 2009Willott et al. , 2010bWillott et al. , 2010aCurtis-Lake et al. 2012;Mallery et al. 2012;Masters et al. 2012;Le Fèvre et al. 2013;Willott et al. 2013;Kriek et al. 2015;Bañados 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;Ono et al. 2018;Pâris et al. 2018;Pentericci et al. 2018;Matsuoka et al. 2018aMatsuoka et al. , 2018bShibuya et al. 2018b;Matsuoka et al. 2019;Zhang et al. 2020;Harikane et al. 2020aHarikane et al. , 2020bEndsley et al. 2021;and Garilli et al. 2021). We adopt their classifications between the galaxies and the 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èvre et al. 2013) and the VIMOS Ultra Deep Survey (Tasca et al. 2017), we take into account sources whose reliabilities of the redshift determinations are >70%-75%, i.e., sources with redshiftreliability 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.
In 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 the 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 e = 0.2 are adopted. We use the Calzetti et al. (2000) dust extinction law with reddening of E(B − V ) = 0.16. The intergalactic medium (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. 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 Sérsic index of n = 1.5, and 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 e = 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  β UV  −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 . 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). Figure 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

Selection Completeness and Redshift Distribution
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 caused by z = 5.7 LAE 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 a lower IGM attenuation than Inoue et al. (2014), resulting in lower redshift dropouts, 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.

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 . 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. relatively shallow. We evaluate contamination fractions in our dropout samples with the following three methods.
The 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  24 mag, because the number density of high-redshift galaxies decreases to brighter magnitudes, while that of the 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 multiband 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 c c < star 2 galaxy 2 (see Laigle et al. 2016). We also check the COSMOS Hubble Advanced Camera for Surveys 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.
The 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 lowor 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 ∼ 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.
The third method is a simulation with shallower data, in the same manner as Ono et al. (2018). We use a shallower data set whose depth is comparable with that of the Wide layer in the UD-COSMOS field, the Wide-layer-depth COSMOS data. We assume that the UD-COSMOS 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 radius of 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.
We 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. (2021, in preparation). 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+HSC-SSP 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:U -dropouts For the BX and BM galaxies, we adopt the following color criteria, respectively: For the ugr filter set, we define selection criteria by comparing the positions of stars in the ugr and u * gr diagrams:U-dropouts Selection 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. (2021, in preparation). 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. (2021, in preparation).

Derivation
We 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, where C(m, z) is the selection completeness, i.e., the probability that a galaxy with an apparent magnitude m at redshift z is detected and satisfies the selection criteria, and dV(z)/dz is the differential comoving volume as a function of redshift. The space number densities of the dropouts that are corrected for incompleteness and contamination effects are obtained by calculating where 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 calculate 1σ uncertainties of the space number densities of dropouts, we consider uncertainties of the surface number densities and the contamination fractions.
We 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¯: where 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 Å. For the apparent magnitude m, we use a magnitude in a band whose central wavelength is the nearest to the rest-frame wavelength of 1500 Å, namely i, z, y, and y bands for g-, r-, i-, and z-dropouts, respectively. We set the Kcorrection 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.

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  −24 (e.g., Akiyama et al. 2018), studies of galaxies at M UV  −22 (e.g., Finkelstein et al. 2015b;Bouwens et al. 2021), and studies of galaxies and AGNs (Ono et al. 2018;Stevans et al. 2018;Adams et al. 2020).
Our derived luminosity functions at z ∼ 5 and 6 are shown in the bottom panel of Figure 5 and the top panel of Figure 6,  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.
respectively. Similar to the z ∼ 4 result, our measurements agree well with previous studies of quasars at M UV  −24 mag (e.g., Matsuoka et al. 2018c;Niida et al. 2020), studies of galaxies at M UV  −22 mag (e.g., Finkelstein et al. 2015b;Bouwens et al. 2021), and studies of galaxies and AGNs (Ono et al. 2018). At z ∼ 7, our derived luminosity function agrees with previous studies (e.g., Bowler et al. 2017;Ono et al. 2018), as shown in the bottom panel of Figure 6. Table 5 summarizes our measurements of the luminosity functions at z ∼ 4-7.
These 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 spectroscopically identified galaxies and AGNs. Based on our spectroscopic results and the literature, we derive the galaxy fraction that 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).

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 Table 5 Obtained Dropout (Galaxy+AGN) and Galaxy UV Luminosity Functions at z ∼ 4, 5, 6, and 7    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: where f * is the overall normalization, L * is the characteristic luminosity, and α and β are the faint and bright-end power-law slopes, respectively. We define a DPL function as a function of where * M UV is the characteristic magnitude. For the galaxy luminosity function, we fit with a DPL function or the Schechter function (Schechter 1976): where f * , * L UV , and α are the overall normalization, the characteristic luminosity, and the faint power-law slope, respectively. We define the Schechter function as a function  Figure 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 f * parameter by 0.7 dex following f * ∝ 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  M UV −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., Akiyama et al. 2018;Ono et al. 2018;Matsuoka et al. 2018c;Niida et al. 2020). Figure 8 summarizes UV luminosity function estimates at z ∼ 4-7 in this work and the literature, and the best-fit 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  −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   This indicates a very rapid growth of AGNs in the first 1.5 Gyr.

Redshift Evolution
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  −26 mag) will be very small (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).

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 : Note 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 compare the estimated galaxy fractions at z ∼ 4 with those derived from the VVDS data (Le Fèvre 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 accurately determine the spectroscopic galaxy fractions. Figure 10 and Table 5 show our estimates of the galaxy UV luminosity functions at z ∼ 4-7. We confirm that our results are consistent with the previous results in the UV magnitude range fainter than −22 mag. This is because the number density of Figure 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, gray, 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.   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).

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., Takahashi et al. 2011;Wyithe et al. 2011;Barone-Nugent et al. 2015, Mason et al. 2015b). 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 stronglensing 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 demagnification of unlensed sources: is 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 ray-tracing 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.
In 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 best-fit 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)σ, and 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. 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 M 21 UV UV 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  M UV  −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 ∼ 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.

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 UV th ). The number of dropouts in the subsamples and their magnitude thresholds are summarized in Table 8. We do not use sources brighter than = m 20.0 UV cut 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.
We calculate observed angular correlation functions of the subsamples, ω obs (θ), using an estimator proposed by Landy & Szalay (1993), where DD(θ), DR(θ), and RR(θ) are the numbers of galaxygalaxy, galaxy-random, and random-random pairs 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 . Corrections for contaminations (e.g., Equation (29) in Harikane et al. 2016) are not applied because the clustering strength of the interlopers is not wellmeasured. 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), andHildebrandt et al. (2009).
Due 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 where ω 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. We estimate the statistical errors of the angular correlation functions using the Jackknife estimator. We divide each subsample into Jackknife samples of about 1000 arcsec 2 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 å w q w q w q w q where N Jack is the total number of the Jackknife samples, ω l is the estimated correlation function from the lth realization, and w 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.

Halo Occupation Distribution (HOD) Model Fitting
We use a halo occupation distribution (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., Ma & Fry 2000;Peacock & Smith 2000;Seljak 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 number densities from the HOD model. Details of the calculations are presented in Harikane et al. (2016).
We 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, where -C i j , 1 is the inverse covariance matrix, n g is a space number density of galaxies in the subsample, and s n log 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 parameter estimation technique. In our HOD model, an occupation function for central galaxies follows a step function with a smooth transition, An occupation function for satellite galaxies is expressed by a power law with a mass cut, The total occupation function is These 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;Garel et al. 2015;Zheng et al. 2005). Indeed, previous studies demonstrate that this HOD model can explain observed angular correlation functions of high-redshift galaxies (Harikane et al. 2016(Harikane et al. , 2018aIshikawa et al. 2017). We calculate the mean dark matter halo mass of galaxies including both the central and satellite galaxies, á ñ M h , effective galaxy bias, b g eff , and the satellite fraction, f sat , as follows:     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(Navarro et al. , 1997, the Duffy et al. (2008) concentration parameter, and the Smith et al. (2003) nonlinear matter power spectrum.
Some theoretical studies claim that the halo bias is scaledependent in the quasi-linear scale of r ∼ 50 Mpc (the nonlinear halo bias effect; Reed et al. 2009;Jose et al. 2013Jose et al. , 2016Jose et al. , 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 correlation functions at 10″ < θ < 90″ (10″ < θ < 120″) in the UltraDeep and Deep (Wide) layers, because they could be affected by the nonlinear halo bias effect. We also do not use the measurements at θ 2″ that are possibly affected by the source confusion.
The HOD model has five parameters, M min , s M log 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 s = 0.2 M log h 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 from M min , we use the relation which 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 We 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  θ  10″) and large (θ  100″) scales. However, the models underpredict the correlation functions by a factor of 1.5-6 in 10″  θ  100″, the transition scale between 1-and 2-halo terms (the quasi-linear scale), especially in the subsamples at z  3. These results indicate that the correlation functions at 10″  θ  100″ cannot be explained by the scale-independent halo bias due to the nonlinear halo bias effect in this quasi-linear scale (Jose et al. 2013(Jose et al. , 2016(Jose et al. , 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 , s M log h , M sat , and α as free parameters, but the results does not significantly change. In Figure 15, we compare the observed number densities and the predictions from the best-fit HOD models, showing good agreement.  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 UV th as M h and M UV , respectively. The crosses are results of the previous work based on the early HSC-SSP data ). The solid curves show the best-fit relations of Equation (54). Figure 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 UV th . Figure 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(Harikane et al. , 2018a at z ∼ 4-7. We plot M min and M UV th 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 e to 1 × 10 13 M e , 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. 2013Behroozi et al. , 2019Moster et al. 2013Moster et al. , 2018Coupon 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. 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 with increasing redshift at the 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 DPL function:  (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.

M UV − M h Relation
In 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 recalculated 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 data sets. 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 previous 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. Note.
(2) Threshold absolute magnitude in the rest-frame UV band. Furthermore, our results improve on both the statistics and the dynamic range covered in M UV .

5.4.
 We estimate a ratio of the SFR to the dark matter accretion rate,  M SFR h , or the baryon conversion efficiency, /Ω 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 , not the ratio of the SFR to the gas mass (SFR/M gas ), which is usually used in radio astronomy. We derive the dust-uncorrected SFRs (SFR UV ) from UV luminosities using the calibration used in Madau & Dickinson (2014) with the Salpeter ( We 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  4 and linearly extrapolate the relation with fixing the slope at z  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 Nbody simulation results in Behroozi & Silk (2015, their 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.
We plot the dust-corrected  M SFR h ratios at z ∼ 2-7 as a function of the halo mass in Figure 18. The results are also summarized in Table 9. The black solid curve in Figure  This relation agrees with the measured  M SFR 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).
On the other hand, the  M SFR 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, 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 recalculated 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 residing in massive halos (see Table 8). ( )) 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  M M SFR h 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.
This 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.
In Figure 19, we compare our  M M SFR h h relation with the results in the literature (Behroozi et al. 2013;Mason et al. 2015a;Moster et al. 2018;Tacchella et al. 2018;Harikane et al. 2018a;Behroozi et al. 2019). The relations in the literature show similar trends to our result; the  M SFR h ratio has a peak of M SFR 0.1 0.01 h -around the halo mass of 10 11 -10 12 M e . However, the slopes of the relation at the high-mass and lowmass ends are different between these studies. These differences are possibly due to differences in used observational data sets, the halo mass functions, and details of modeling. We also note that there are several systematic uncertainties on the  M SFR 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 (Bouwens et al. , 2020Finkelstein et al. 2015b;Oesch et al. 2018).
First we assume the redshift-independent  M M SFR h h relation (Equation (56)).   (55)). The result of Harikane et al. (2018a) is recalculated based on the accretion rate of Behroozi & Silk (2015), and thus is ∼2 times higher than the original relation in Harikane et al. (2018a). 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.
Then we use the gradually evolving  M M SFR h 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  M M SFR h 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  M M SFR h h relations by using normal star-forming galaxies (dropout, BX, and BM galaxies), and not considered quiescent or dusty starburst galaxies. Quiescent and dusty starburst galaxies are expected to have lower and higher  M SFR h ratios than the normal star-forming galaxies, and effects of these galaxies are thought to be less dominant at z  4 (e.g., Muzzin et al. 2013;Bouwens et al. 2020;Fudamoto et al. 2021). Although quiescent and dusty starburst galaxies would be nonnegligible at z  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.
The 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  M M SFR h 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  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  4. The dark matter (and gas) accretion rate monotonically decreases over the whole redshift range 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 certain 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   -M M SFR 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.
In our calculation, we integrate Equation (58) down to the SFR of 0.3 M e yr −1 , corresponding to the halo mass of ∼3 × 10 10 M e at z ∼ 7. This integration limit is chosen to match the calculations at z  2 in previous studies (Bouwens et al. , 2020Finkelstein et al. 2015b;Oesch et al. 2018), and slightly different from the calculations in Madau & Dickinson (2014), who integrate down to 0.03L * . This difference does not affect our discussions above. Bouwens et al. (2021) show that the * M UV parameter of the luminosity Figure 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 e , 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).  (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.  (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. 2016Harikane et al. , 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 the observed results given the 0.15 dex (1σ) uncertainty in  M SFR h , indicating that our  M M SFR h h relation is consistent with the observed redshift evolution of the UV luminosity function. Mason 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  4. This is consistent with our results of the  M M SFR h 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 density, as discussed above.

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  z  10. By assuming that this  M M SFR h 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 the dark matter accretion rate. Figure 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;Oesch et al. 2018;Tacchella 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  10 and our predictions at z  10 are well fitted with the following function:  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). The 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. Using this method to predict the SFR density at z > 10, we can obtain a rough estimate of the epoch of the first star formation. We calculate the cumulative number of formed stars as a function of redshift: where 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 e (e.g., Hirano et al. 2015). We adopt = -V h 3 Mpc survey 1 3 ( ) that is the volume of the simulation box in Hirano et al. (2015). We calculate ρ SFR using Equation (58) by extrapolating the  M M SFR h h relation both to the higher-redshift and lowermass range. We integrate down to the halo mass of 10 5 M e , comparable to the halo masses of the first stars in the 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  M M SFR h h relation and the mass limit of the integration that we adopt 10 4 -10 6 M e . 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(Hirano et al. , 2015, although this is a very rough estimate. In particular, it is not clear whether the assumed  M M SFR h 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.

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  −23 mag), indicating the existence of the bright-end excess of the number density beyond the Schechter function. Since these luminosity Figure 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 e yr −1 (M UV = −17 mag), as previous studies (Bouwens et al. , 2020Finkelstein et al. 2015b;Oesch et al. 2018). The gray dashed curve shows the extrapolation of the relation of Madau & Dickinson (2014) (60)) to the observed cosmic SFR densities at z  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)).
Then we calculate the SFR and mass accretion rate in the same manner as Section 5.4, and obtain the  M SFR h ratio. Figure 26 presents the luminosity function, the M UV -M h relation, and the star formation efficiency (the  M M SFR h 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 (  M SFR 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 ; 10 13 (10 12 ) M e 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 timescale 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. (D) Low dust obscuration. Bowler et al. (2020) discuss 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;Bakx et al. 2020;Fudamoto et al. 2020;Harikane et al. 2020b), and the β UV -M UV relation is not wellconstrained at this very bright magnitude range (i.e., ∼−23 mag). In addition, some recent studies report dustpoor UV-luminous star-forming galaxies at z > 2 (Marques-Chaves 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. (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  −22 mag) galaxies at z  7 (e.g., Laporte et al. 2017;Mainali et al. 2018;Endsley et al. 2021;Jiang 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.
Based 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 redshift. We cannot distinguish these possibilities with the current data sets, and future large and deep observations are needed. For example, the 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.

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.
Our major findings are summarized below: 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  M UV −14 mag, corresponding to the luminosity range of  * L 0.002 UV  * L L 2000 UV UV . 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 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). 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 2σ levels ( Figure 10), which is possibly made by the inefficient mass quenching, low dust obscuration, and/or hidden AGN activity (Section 6.3). 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). 5. We have calculated the ratio of the SFR to the dark matter accretion rate,  M SFR h , and identified an  M M SFR h h relation that 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  M SFR h ratio gradually increases with decreasing redshift from z ∼ 5 to 2 within 0.3 dex.