Black Hole Mass and Eddington-ratio Distributions of Less-luminous Quasars at z ∼ 4 in the Subaru Hyper Suprime-Cam Wide Field

We investigate the black hole mass function (BHMF) and Eddington-ratio distribution function (ERDF) of broad-line active galactic nuclei (AGNs) at z = 4, based on a sample of 52 quasars with i < 23.2 at 3.50 ≤ z ≤ 4.25 from the Hyper Suprime-Cam Subaru Strategic Program S16A-Wide2 data set, and 1462 quasars with i < 20.2 in the same redshift range from the Sloan Digital Sky Survey data release 7 quasar catalog. Virial black hole (BH) masses of quasars are estimated using the width of the CIV 1549 Å line and the continuum luminosity at 1350 Å. To obtain the intrinsic broad-line AGN BHMF and ERDF, we correct for the incompleteness in the low-mass and/or low-Eddington-ratio ranges caused by the flux-limited selection. The resulting BHMF is constrained down to logMBH/M⊙∼7.5 . In comparison with broad-line AGN BHMFs at z ∼ 2 in the literature, we find that the number density of massive SMBHs peaks at higher redshifts, consistent with the downsizing evolutionary scenario. Additionally, the resulting ERDF shows a negative dependence on BH mass, suggesting more massive SMBHs tend to accrete at lower-Eddington ratios at z = 4. With the derived intrinsic broad-line AGN BHMF, we also evaluate the active fraction of broad-line AGNs among the entire SMBH population at z = 4. The resulting active fraction may suggest a positive dependence on BH mass. Finally, we examine the time evolution of broad-line AGN BHMF between z = 4 and 6 through solving the continuity equation. The results suggest that the broad-line AGN BHMFs at z = 4–6 only show evolution in their normalization, but with no significant changes in their shape.


INTRODUCTION
Understanding the formation and evolution of supermassive black holes (SMBHs), which are found to be ubiquitous at the centers of local massive galaxies (e.g., Kormendy & Richstone 1995), is one of the important goals in observational cosmology.The growth processes of SMBHs may be physically associated with the galaxy evolution processes, as suggested by a series of tight scaling relationships observed between the mass of SMBHs and the properties of the host spheroids, such as their stellar velocity dispersion (e.g., Gebhardt et al. 2000;Merritt & Ferrarese 2001), bulge mass (e.g., Magorrian et al. 1998;Marconi & Hunt 2003).Those empirical correlations imply a physical connection between the evolution of SMBHs and host galaxies, i.e., the "coevolution".But its physical mechanism is still under debate: in one scenario, galaxy minor mergers are proposed to naturally invoke the hierarchical assembly of SMBHs and stellar mass of host galaxies (e.g., Peng et al. 2007;Jahnke et al. 2011); in another scenario, galaxy major mergers and the following AGN feedback are thought to be the mechanism to regulate the growth of SMBHs together with star formation in the central region of their host galaxies (e.g., Fabian 1999;Di Matteo et al. 2005;Springel et al. 2005;Hopkins et al. 2006).
Luminosity functions of AGNs can reflect the SMBH growth under the active accretion phase, and they are already intensively investigated using large samples of AGNs across cosmic time (e.g., Ueda et al. 2003;Fan 2004;Richards et al. 2006a;Silverman et al. 2008;Croom et al. 2009;Ross et al. 2013;Toba et al. 2013Toba et al. , 2014;;Ueda et al. 2014;Akiyama et al. 2018;McGreer et al. 2018;Matsuoka et al. 2018Matsuoka et al. , 2023)).Following the idea proposed by Soltan (1982), the total mass accreted towards SMBHs, which is obtained by integrating the AGN luminosity function through luminosity and redshift, is compared to the mass density of SMBHs in the local universe (e.g., Yu & Tremaine 2002;Shankar et al. 2009).Consistency between the two quantities indicates the local SMBH mass density is mostly accumulated by mass accretion in the AGN phase, and requires the average radiative efficiency to be in the order of 10%.
Based on luminosity functions, the cosmological evolution of the comoving number density of AGNs has been examined.The number density of luminous AGNs, e.g., quasars, is found to peak at higher redshifts than that of less-luminous AGNs (e.g., Ueda et al. 2014;Akiyama et al. 2018;Matsuoka et al. 2018;Niida et al. 2020).Such anti-hierarchical trend is called "down-sizing" evolution.At high redshifts, there are also observational results suggesting that the abundance of less-luminous AGNs tends to decline more mildly towards higher red-shifts than that of luminous AGNs (e.g., Ueda et al. 2014).Such trend is called "up-sizing" evolution, which follows the direct expectation of hierarchical structure formation.
However, the luminosity of AGNs does not directly reflect the mass assembly history of SMBHs.The luminosity function of AGNs is the convolution of the black hole mass function (BHMF) and the Eddington-ratio distribution function (ERDF).Some studies find that the Eddington ratio, a ratio between the observed bolometric and the Eddington-limited luminosity, is distributed over a wide range (e.g., Schulze & Wisotzki 2010;Suh et al. 2015;Jones et al. 2016).Thus, constraining BHMF and ERDF along the cosmic time becomes important in studying the growth history of SMBHs quantitatively.In the local universe, by applying the scaling relationships between the mass of SMBHs and properties of their host spheroids for normal galaxies, the BHMF is evaluated in many studies (e.g., Yu & Tremaine 2002;Marconi et al. 2004;Shankar et al. 2009).Beyond the local universe, since it remains controversial how the scaling relationships evolve in the early universe, determining the BHMF by using the scaling relationships is not applicable.
With the BH mass estimates, some studies jointly determine the BHMF and ERDF by fitting bivariate distribution models of BH mass and Eddington ratio with observed distributions, either using the Bayesian method (e.g., Shen & Kelly 2012;Kelly & Shen 2013), or the Maximum-likelihood method (e.g., Schulze & Wisotzki 2010;Nobuta et al. 2012;Schulze et al. 2015;Ananna et al. 2022).While different methods and quasar samples are used, those studies result in a common trend of "down-sizing" in the growth history of SMBHs: the number density of quasars with less-massive SMBHs mildly increases from z = 2 − 3 to the local universe, while that of quasars with massive SMBHs sharply declines towards z = 0, i.e., the massive SMBHs stop growing first while the less-massive ones grow continuously.
The evolutionary trend remains uncertain towards higher redshifts.Most of the current statistical studies only rely on the luminous quasars brighter than i = 20 at z ≳ 3. Due to the flux limit, a large fraction of AGNs with low Eddington ratios of log λ Edd < −0.5 or small BH masses of log M BH /M ⊙ < 9 can be missed, and the constraints on the less-massive end of BHMF and low-Eddington-ratio end of ERDF are not tight.At z > 3.2, it is found that the fraction of observable quasars above the SDSS detection limit decreases to less than 10% at M BH < 10 9 M ⊙ , indicating the constraints can have large uncertainty in the less-massive ends (Shen & Kelly 2012).Schulze et al. (2015) also demonstrate a sample of AGNs covering a wide luminosity range is crucial to constrain the BHMFs over a wide BH mass range.
In this work, in order to trace the early growth of SMBHs, we constrain the BHMF and ERDF of broadline AGNs at z = 4 using a sample of less-luminous quasars selected from the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; Aihara et al. 2018a).Here, HSC-SSP is an imaging survey utilizing Hyper Suprime-Cam (HSC), the wide-field CCD camera attached to the prime focus of the Subaru 8.2m telescope (Miyazaki et al. 2018).In the most recent data release PDR3, the Wide layer of HSC-SSP covers ∼1200 deg 2 in the grizy-bands with the 5σ detection limits of 26.5, 26.5, 26.2, 25.2, and 24.4 mag, respectively (Aihara et al. 2022).Using the wide coverage and deep detection limit of the HSC-SSP, Akiyama et al. (2018) construct a sample of 1666 z ∼ 4 less-luminous quasars from the S16A-Wide2 data release, which has a sky coverage of ∼ 170 deg 2 (Aihara et al. 2018b).The sample covers ∼ 2 − 3 mag fainter magnitude range than the SDSS quasars at z = 4, with the sample size a few times larger than other z ∼ 3 − 4 quasar samples in the similar luminosity ranges (e.g., Netzer et al. 2007a;Trakhtenbrot et al. 2016;Schulze et al. 2018).
The outline of this paper is as follows.We describe the quasar samples and spectroscopic observations in section 2. Line fitting of the C IV broad emission line, and estimates of the virial BH mass and Eddington ratio are described in section 3. Section 4 introduces the V max and Maximum-likelihood method used to constrain the z = 4 broad-line AGN BHMF and ERDF.The final section 5 discusses the cosmological evolution of the broad-line AGN BHMFs and ERDFs.In addition, we investigate the active fraction of broad-line AGNs at z = 4 based on the derived BHMFs and ERDFs.We also examine the evolution of broad-line AGN BHMFs and ERDFs between z = 4 and 6 through solving the continuity equation.Throughout this paper we assume a ΛCDM cosmology with Ω m = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 .The PSF magnitudes of quasars, which are determined by fitting a model PSF to the image of object, are adopted.All the magnitudes are given in the AB magnitude system.We use the sample of less-luminous z ∼ 4 quasar candidates constructed with the HSC-SSP S16A-Wide2 dataset (for details see Akiyama et al. 2018).To build the z ∼ 4 quasar catalog, Akiyama et al. (2018) adopt the following g-band dropout color criteria: 0.65(g − r) − 0.30 > (r − z), (1) 3.50(g − r) − 2.90 > (r − z), (g − r) < 1.5. (3) Additional color criteria are set to remove contamination from Galactic stars as follows: (i − z) > −0.3. (5) We refer to the color criteria as "c1".By applying the c1 criteria for objects with stellar morphology (for the criteria see equation 1-2 in Akiyama et al. 2018) on the HSC-SSP S16A-Wide2 dataset, a sample of 1,666 z ∼ 4 quasar candidates at 20 < i < 24 is constructed in an effective area of 172 deg 2 .Their photometric redshifts are estimated to distribute between z phot = 3.4 and 4.6, with a mean and standard deviation at z phot = 3.9 and 0.2, respectively (Akiyama et al. 2018).
Additionally, in order to examine the completeness of the c1 criteria for selecting quasars at 3.5 < z < 4.5, we include supplemental z ∼ 4 quasar candidates in the spectroscopic observations.Since the c1 criteria are determined to minimize the contamination from compact galaxies, Galactic stars, and quasars at other redshifts than 3.5 < z < 4.5, they can suffer from a relatively low completeness.Thus, we expand the color selection window on the g − r vs. r − z plane to recover the z ∼ 4 less-luminous quasars that are missed by the c1 criteria.
The additional c2-5 selection criteria are summarized in Table 1.In Figure 1, the red, blue, magenta, yellow and green shaded areas show the c1-5 selection criteria,  2018), with black squares implying its colors at z = 2, 2.5, 3, 3.5, 4, and 4.5 from left to right.Stars and dots show the identified quasars with QOP=4 (good quality) and QOP=3 (low quality), respectively.Symbols are color-coded with their spectroscopic redshift, as suggested by the color bar on the right side.Black open circles mark quasars selected by the Lyman break criteria in Ono et al. (2018) or radio detection in Yamashita et al. (2018).Right) the g − r vs. r − z color-color diagram of the identified contaminants.Orange squares, green triangles and purple inverted triangles represent the z < 3 quasars, galaxies and Galactic stars identified in the spectroscopic observations, respectively.
Table 1.The c1-5 color selection criteria and number of the observed candidates in the spectroscopic observations respectively.Here, the selection windows are designed according to the mean color track of the model quasars in Akiyama et al. (2018), which is plotted as a function of redshift by solid black line.Since the model quasars do not consider the dust reddening to the nucleus, we set the c2 criteria to include quasars at the same redshifts with c1 but reddened by dust.The region also contains quasars in the redshift range of 4.0 < z < 4.5, where the completeness of the c1 criteria significantly drops.The c3 criteria extend to the area where con-taminations by quasars at z < 3.5 and Galactic stars dramatically increase.The c4 and c5 criteria further expand the selection area to cover those z ∼ 4 quasars with large offsets away from the model colors.By applying the supplemental c2-5 color criteria for objects with stellar morphology in the HSC-SSP S16A-Wide2 dataset, we select additional 11009 c2, 10060 c3, 84412 c4 and 51525 c5 targets with 20 < i < 24.They are included as fillers for the multi-object spectroscopic observations.

Spectroscopic observations and data processing
We carried out spectroscopic observations of the z ∼ 4 quasar candidates using the two-degree field (2dF) fibre positioner (Lews et al. 2002) with the AAOmega spectrograph (Smith et al. 2004) mounted on the 3.9m Anglo-Australian Telescope (AAT) during three observing runs S17A/08, S18A/16, and S18B/03 between 2017 to 2019, and the DEep Imaging Multi-Object Spectrograph (DEIMOS) (Faber et al. 2003) mounted on the Keck II telescope during the observing run S356 in 2019.

Observations with AAT/AAOmega
The 2dF fibre positioner is able to place up to 400 fibres (8 fibre bundles for guiding, 25 fibres for sky subtraction, and remaining fibres for science objects) within a 2 deg diameter FoV, and AAOmega is a dual-beam spectrograph.In the AAT S17A/08 observing run, to fully cover both of the Lyα 1215 Å and C IV 1549 Å emission lines of quasars at 3 < z < 4.5, we choose the follow-ing setup: the X6700 dichroic beam splitter, the 580V grating with central wavelength at 5750 Å in the blue channel, and the 385R grating with central wavelength at 8200 Å in the red channel.The configuration covers a wavelength range of 4800-9800 Å, with a resolution of R∼1600.Target FoVs are determined to maximize the number of c1 quasar candidates, and c2-5 targets are added as fillers in each field.In total, we pick up 7 target fields, which contain 84 c1 and 1695 c2-5 candidates with 20 < i psf < 23.During the observing run, the seeing is ∼ 1.4 ′′ − 1.8 ′′ .Each field is observed with a net exposure time of 160-300 minutes, typically divided into multiple exposures of 20-30 minutes.Details are summarized in Table 2.
In the AAT S18A/16 and S18B/03 observing runs, since a spectral break between the red and blue channels after the splicing process is found in some spectra with the setup using the dichroic beam splitter at 6700 Å, we choose the standard setups: the X5700 dichroic beam splitter, the 580V grating with central wavelength at 4821 Å in the blue channel, and the 385R grating with central wavelength at 7251 Å in the red channel.The configuration covers a wavelength range of 3800-8800 Å, with a resolution of R∼1400.In the S18A/16 run, we observe one FoV containing 16 c1 and 236 c2-5 candidates with 20 < i psf < 23.5.In the S18B/03 observing run, we observe three FoVs including 93 c1 and 52 c3 candidates with 20 < i psf < 24 as fillers of a program targeting z > 4 Lyman break galaxies.Both of the observing runs have bad weather conditions, and the effective exposure time is severely reduced, resulting in low SNRs in most of the spectra.Details of the S18A/16 and S18B/03 observing runs are summarized in Table 2 and 3, respectively.
The raw data is processed by the OzDES pipeline, which uses python scripts to reduce the individual science frame with the v6.46 2dfdr pipeline, and then to combine the spectra from multiple exposures (private communication with Lidman, C.).Wavelength dependence of the sensitivity is corrected using the default sensitivity function provided in the pipeline.PyCosmic (Husemann et al. 2012) is used to remove pixels affected by cosmic rays.There are 5 spectra affected by a severe splicing break at 6700 Å between the red and blue channels, thus we manually re-normalize the blue channel by matching the mean count in 6600-6700 Å in the blue channel to that in 6800-6900 Å in the red channel.For each spectrum, flux calibration is applied by adjusting its normalization to minimize the average difference between the riz-band PSF magnitudes cataloged in the HSC-SSP S16A-Wide2 dataset, and those directly calculated from the reduced spectrum.The riz-bands are set since they are mostly covered by the observed wavelength range.The cataloged magnitudes have been corrected for the galactic extinction.

Observation with Keck/DEIMOS
DEIMOS is a multi-slit imaging spectrograph, which covers up to ∼ 100 objects within a FoV of 16 ′ × 4 ′ .In the S19A/S356 observing run, we adopt a setup with a slit width of 1 ′′ , the 600 line mm −1 (600ZD) grating with central wavelength at 7500 Å, and the GG495 blue blocking filter to cover both of the Lyα 1215 Å and C IV 1549 Å emission lines at 3 < z < 4.5.The configuration covers a typical wavelength range of ∼5000-10000 Å with a resolution of R ∼1600, but the range depends on the target position on the mask.We observe 12 FoVs, which contain in total 41 c1 and 144 c2-5 candidates with 20 < i psf < 23.5.The seeing condition is ∼ 0.5 ′′ − 0.7 ′′ .Each field is observed with a total exposure time of 40-90 minutes, typically divided into 2-3 exposures of 20-30 minutes.Details are summarized in Table 2.
The raw data is processed by the DEEP2 pipeline (Cooper et al. 2012;Newman et al. 2013).At first bias subtraction, cosmic ray removal, and flat-fielding are applied in a standard manner, then wavelength calibration and slit-tilt correction are determined.Sky subtraction is applied to the tilt-corrected data, and extractions of spectra are conducted.The extracted spectra from multiple frames are combined.Two CCD chips 2 and 6 of the J091408−012722 field, covering one c1 and one c4 candidates, fail in the reduction with the pipeline.We then use the IRAF longslit and twodspec packages to reduce their spectra by cutting the data into longslit spectra.Feige 34 and BD+25d3941 are observed to obtain the throughput as a function of wavelength, and the flux calibration is applied using the same method described in the AAOmega data processing.
In total, we cover 234 c1, 216 c2, 324 c3, 794 c4, and 793 c5 candidates, as summarized in Table 1.We observe 141 c1 quasars on the HSC Wide layer, which account for 141/1666 = 8.5% of the full quasar sample in Akiyama et al. (2018).The photometric redshift of the 141 quasars spans between z phot = 3.4 and 4.5, with the mean at z = 3.8 and the standard deviation of 0.2.This distribution is consistent with that of the full quasar sample in Akiyama et al. (2018), suggesting we do not introduce bias in the distribution of photometric redshift when selecting targets for spectroscopic observations.

Classification and redshift determination
We use the web-based program MARZ (Hinton et al. 2016) to manually inspect each spectrum and determine its redshift.MARZ can automatically determine the bestfit object type and redshift.It matches input spectra against a variety of spectral templates of stars and galaxies using a modified version of AUTOZ cross-correlation algorithm implemented by Baldry et al. (2014).For an incorrect matching, users are allowed to manually re-fit the spectrum by iterating the automatic results, modifying the choice of a template, or marking spe-cific spectral features.MARZ is originally developed for the Australian Dark Energy Survey (OzDES) carried out with the AAT/2dF-AAOmega spectrograph, while it can also handle spectrum taken with other instruments and pipelines if the spectrum is in the standard FITS format.In the course of inspection, we assign a quality flag to each spectrum with the following criteria: • QOP=4 -quasars at 3 < z < 4.5 with both the Lyα and C IV emission lines detected, and the C IV line profile detected with sufficient SNRs; • QOP=3 -quasars at 3 < z < 4.5 with only the Lyα emission line detected, but the C IV emission line highly affected by the absorption lines, or barely detected that its profile can not be well constrained; • QOP=2 -quasars at z < 3 with other multiple emission lines detected; • QOP=1 -galaxies; • QOP=6 -stars; • QOP=0 -unknown or non-detection.
Two example spectra assigned with the quality flag of QOP=4 and 3 are shown in the top and bottom panels of Figure 2, respectively.There are 125 out of the 234 c1 candidates observed under good observing conditions, i.e., during the S17A/08 and S19A/S356 runs.Among them, we identify 78 quasars at 3 < z < 4.5 with 51 QOP=4 and 27 QOP=3 quality flags.The resulting success rate of c1 candidates to be a z ∼ 4 quasar is 78/125 = 62.4%.There are 1 and 4 c1 candidates identified as a Galactic star and z < 3 quasars, respectively.The remaining objects show no clear features to be identified with any of the templates due to either low SNRs or failure in the reduction of their spectra.Among the 150 c2, 212 c3, 744 c4 and 733 c5 candidates observed in the two observing runs, 9, 18, 6 and 3 of them are identified with a 3 < z < 4.5 quasar, resulting in the identification rates of 6%, 8.5%, 0.8% and 0.4%, respectively.
In total, we identify 151 quasars at 3 < z < 4.5, containing 80 with QOP=4, and 71 with QOP=3 quality flags.Among them, 11 show strong and broad absorption feature of the C IV line, and they are identified as broad absorption line (BAL) quasars.Details are summarized in Table 5.
In Figure 1, we plot the g − r and r − z color distribution of the identified sources.The identified 3 < z < 4.5 quasars are shown by stars (QOP=4) and dots (QOP=3) in the left panel.Their distribution follows the color track of model quasars between z = 3 and 4.5.There are 9 quasars, marked by open black circles, locating outside of the c1-5 windows due to different selection criteria of the Lyman break color or radio detection.Additional 5 quasars fall in the c1-3 selection windows, but they are removed from our samples by the i−z and z −y color criteria.In the right panel, color distributions of the contaminating z < 3 quasars, galaxies and Galactic stars are plotted by orange squares, green triangles and purple inverted triangles, respectively.Majority of them locate in the c3-5 selection windows as expected, with a few contaminants locating at the edges of the c1 window.The 4000 Å break in z ∼ 0.5 galaxies can mimic the Lyman break at z ∼ 4. Also, the broad Mg II 2798 Å emission line in some quasar spectra mimics the Lyα 1215 Å emission line.
In Figure 3, we plot the i-band magnitude vs. redshift distribution of the identified quasars at 3 < z < 4.5.The pink open histogram in the horizontal panel shows that the 151 identified quasars cover a redshift range of 3.04 ≤ z spec ≤ 4.44, with a mean and standard deviation at 3.67 and 0.28, respectively.The vertical panel displays the i-band magnitude distribution of the identified quasars.Pink open histogram represents the entire 151 identified quasars, and they cover a magnitude range of 18.9 < i psf < 23.2, with a median at i = 21.7.If only considering the QOP=4 quasars, whose distribution is plotted by red open histogram, we see they cover a similar magnitude range, but with a slightly brighter median at i = 21.5, indicating no or weak detection of the C IV emission line in the QOP=3 quasars is partly due to their faintness.
Here, we examine the completeness of the c1 criteria in selecting quasars at 3.5 ≤ z ≤ 4.5.Since the c1-5 criteria are determined following the two-color distributions of model quasars at 3.5 < z < 4.5, we assume the combined c1-5 criteria are wide enough to cover the typical quasars in the redshift range.We use the data taken in the observing run S17A/08 as it covers a large number of c1-5 candidates under good conditions.There are 84 c1, 138 c2, 198 c3, 682 c4 and 677 c5 candidates observed in total, which occupy 13.2%, 25.0%, 9.4%, 1.3% and 2.0% of the full c1-5 samples down to i = 23, respectively.Among the observed candidates, we identify 60 c1, 9 c2, 7 c3 and 1 c5 candidates to be quasars at 3.5 ≤ z ≤ 4.5.As discussed in section 2.2, we do not introduce bias in photometric redshift distribution when selecting the c1 candidates for spectroscopic observations, and c2-5 candidates are randomly selected over the target fields as fillers.Based on the statistics, the completeness a Numbers in the brackets are for the BAL quasar identified in each criteria.
of c1 criteria in selecting quasars at 3.5 ≤ z ≤ 4.5 is (60/0.132)/(60/0.132+9/0.250+7/0.094+1/0.020)=74%,which is broadly consistent with that evaluated with the SDSS z ∼ 4 quasars (66%; for details see Akiyama et al. 2018).Except for the c1 criteria, the c3 criteria show the highest completeness of 12%.However, the c3 criteria also show much higher contamination rate of 9% than that of the c1 criteria of 1% in the magnitude range of 20 < i < 23.We thus confirm the c1 criteria established in Akiyama et al. (2018) are able to select quasars at 3.5 ≤ z ≤ 4.5 with high completeness and the minimal contamination.
The photometric redshift of c1 quasar candidates are estimated in Akiyama et al. (2018).We iden-tify 93 among them in the spectroscopic observations.The comparison between their photometric and spectroscopic redshifts is shown in Figure 4.The identified quasars with QOP=4 and QOP=3 quality flags are plotted by red stars and pink dots, respectively.Overall, their photometric redshifts are consistent with the spectroscopic redshifts.The accuracy of photometric redshift, which is defined as 1.48 × median(|z phot − z spec |/(1 + z spec )) (e.g., Ilbert et al. 2006), is 0.04.No quasars fall in the outlier of |z phot − z spec |/(1 + z spec ) > 0.15.We note that there are five quasars with the discrepancy between their photometric and spectroscopic redshifts larger than 0.35.After inspecting their spectra, we find the one in the lower right corner is a BAL  Note-Table 5 is published in its entirety in the machinereadable format.A portion is shown here for guidance regarding its form and content.
a Quasars selected by the Lyman break color or radio detection quasar, and the remaining 4 quasars in the upper left corner show flatter continuum shape in their spectra than those expected from their broad-band photometries.The systematic offsets between the flux-calibrated and cataloged magnitudes of the 4 quasars reach 0.5-2 mag, resulting in a large discrepancy between the photometric and spectroscopic redshift.

He et al.
Moreover, we examine the distribution of spectroscopic redshift of the identified 3 < z < 4.5 quasars in each FoV.In the 20 FoVs within the Wide2 dataset, the standard deviation of redshift distribution of each FoV is ∆z spec ∼0.1-0.3, which corresponds to a velocity offset of ∆v=6600-20000 km s −1 .Thus, most of the identified quasars in the target fields are just clustered in the projected plane by chance.There are 2 quasar pairs, which have ∆v<2000 km s −1 and the angular projected distance r ⊥ < 4 proper Mpc (pMpc).They are listed in Table 6.Following the definition of a quasar pair in Onoue et al. (2018), i.e., ∆v<3000 km s −1 and r ⊥ < 4 pMpc, the 2 quasar pairs can have physical association, but none of them meet the tighter criterion for physical association in Chiang et al. (2013), which is determined by the size of simulated proto-clusters at z ∼ 4.

Single-epoch virial BH mass estimator
Under the assumption that the broad-line region (BLR) is virialized, the virial BH mass can be estimated from the motion of the BLR through M BH = f r BLR v 2 /G, where f is a scaling factor accounting for the kinematic structure of the BLR, v is the velocity of the BLR gas that can be derived from line width of broad emission lines, and r BLR is the BLR radius.Based on the reverberation mapping of local broad-line AGNs, a tight correlation between the size of broad Hβ-emitting region (r BLR,Hβ ) and the continuum luminosity around the Hβ emission line at 5100 Å is found (the r − L relationship; Kaspi et al. 2000).The relationship enables us to estimate the virial BH mass of broad line AGNs simply with the Hβ line width (∆v Hβ ) and the 5100 Å continuum luminosity (L 5100 ), both of which can be derived from a single-epoch spectrum.
The virial BH mass estimates have been extended to other broad emission lines, such as Hα, Mg II and C IV emission lines (e.g., Vestergaard & Peterson 2006;Shen & Liu 2012;Denney et al. 2013;Park et al. 2013Park et al. , 2017)).In this work, since only the C IV 1549 Å emission line of the z ∼ 4 quasars is covered among the calibrated broad emission lines in the optical wavelength range, we adopt the C IV virial BH mass estimate in Vestergaard & Peterson (2006) where FWHM represents the line width of C IV 1549 Å emission line, and L 1350 is the monochromic UV luminosity at 1350 Å.The systematic uncertainty of the calibration is evaluated to be 0.36 dex (Vestergaard & Peterson 2006).
We note that there are also other studies providing calibrations of the C IV -based single-epoch virial BH mass estimate (e.g., Shen & Liu 2012;Park et al. 2013Park et al. , 2017)).For example, Park et al. (2013) use a sample of 26 local AGNs with both the Hβ reverberation masses and UV archival spectra available to re-calibrate the C IV -based virial BH mass estimator.Compared to the method used in Vestergaard & Peterson (2006), Park et al. (2013) perform the multi-component fitting on the C IV complex region to precisely deblend C IV from other contaminating lines, and they relax the constraint of the slope parameter in the calibration.The virial BH masses estimated with the calibration suggest those estimated with the calibration in Vestergaard & Peterson (2006) (i.e., equation 6) can be underestimated in the low-mass range of log M BH /M ⊙ < 8 and overestimated in the high-mass range of log M BH /M ⊙ ∼ 9. Park et al. (2017) extend the work by including additional six reverberation-mapped AGNs with BH masses down to 10 6.5 M ⊙ .The calibration remains consistent with the one in Park et al. (2013).Here, since the calibration in Vestergaard & Peterson (2006) is commonly used in literature, we firstly adopt it to estimate the virial BH mass for comparison with literature, and we discuss how the other calibrations can affect the estimates in section 5.5.

Line width and 1350 Å monochromic luminosity
We measure the continuum flux and emission line width for the 80 QOP=4 quasars, whose C IV 1549 Å emission line has sufficient SNRs.Firstly, we fit the continuum of the obtained spectra by a single power-law model F λ ∝ λ α .The fitting is carried out in the rest frame wavelength ranges of 1340-1350 Å, 1450-1460 Å, 1700-1710 Å, 1800-1810 Å and 2000-2010 Å, since the contribution from emission and/or absorption lines in these wavelength ranges is small.We manually inspect the best-fit model of each spectrum, and shift the wavelength ranges for fitting by 10-100 Å to improve fitting results, especially when large systematic residuals in the wavelength range of 1300-1600 Å are found.In the top panel of Figure 2, we show the best-fit power-law model of the continuum of one example spectrum by a red line.
After subtracting the power-law component, the C IV emission line is fitted with Gaussian models over the wavelength range of 1500-1600 Å.We start from a single Gaussian model, and add an additional Gaus-  are 13 C IV line profiles affected by narrow absorption lines.Since the SNRs of the spectra are not high enough to fit both the emission and absorption lines simultaneously, we fit the spectra by excluding the wavelength ranges affected by those absorption lines.In addition, we find large systematic residuals from the continuum subtraction at the edges of the C IV emission line in 15 spectra, and we exclude those wavelength ranges in fitting as well.Moreover, there is one C IV line profile dominated by spiky narrow component.To prevent the influence of the narrow-line region on estimating the virial BH mass, we remove the component from the best-fit model.There are 6 C IV emission lines well fitted with a single Gaussian, but for the remaining 74 spectra, two Gaussians are required to smoothly describe the observed profile.An example of the fitting is shown in Figure 5.
With the best-fit profile constructed by combining all the Gaussian components, we measure the C IV line width.Both the FWHM and line dispersion can be used to describe the broad line width.While there are studies suggesting the line dispersion is more appropriate in measuring the broad line width (e.g., Denney et al. 2013;Park et al. 2017), it requires high S/N data to accurately fit the line wings, especially for the line profile with extended wings like C IV line (e.g., Park et al. 2017).Here, since most of our spectra have limited SNR∼2-5 at the peak of C IV line, and systematic residuals of the background and continuum subtraction are found in some spectra, especially in those taken with AAT/AAOmega, we only measure the FWHM as line width for virial BH mass estimates.The derived FWHMs are corrected for the instrumental broadening of the spectrograph, which is evaluated using the line width of the night sky emission lines.Table 7 summarizes the results.
The monochromic luminosity at 1350 Å, L 1350 , is derived from the best-fit power-law model.Considering the stellar morphology of the identified quasars, we assume contribution of the host galaxy emission is negligible in the UV wavelength range.As described in section 2.2, the spectra taken with the AAT/AAOmega can be affected by a splicing break between the red and blue arms.Thus, we also measure L 1350 using the grizy 5-band PSF magnitudes in the HSC-SSP S16A-Wide2 dataset.The magnitudes have been corrected for the galactic extinction.Firstly we subtract the contribution of the emission lines, such as, Lyα and C IV , to the corresponding broad-band magnitudes with the line flux measured in the continuum-subtracted spectra.Then, we fit the subtracted 5-band magnitudes with a powerlaw model F λ ∝ λ α , and derive L 1350 with the best-fit power-law model.
In Figure 6, we compare L 1350 of quasars derived with their spectra with those derived from broad-band photometries.For the majority of the objects, the two measurements show consistency, with a median discrepancy of −0.01 dex and a standard deviation of 0.15 dex.Therefore, we directly use their L 1350−spec as L 1350 .We  note that there are 7 objects having the discrepancy greater than 0.3 dex.By inspecting their spectra, we find that the two objects below the lower dashed line show weak detection of continuum and are affected by a splicing break.We thus adopt their L 1350−phot as L 1350 afterwards.For the five objects above the upper dashed line, their spectra show large systematic residuals of background subtraction in the wavelength range of λ obs > 6000 Å, resulting in over-subtraction in the izyband photometries.We keep their L 1350−spec as L 1350 .
For the line width, the quasars cover a range of 952 − 10902 km s −1 in FWHM, with a median of 3081 km s −1 .Compared with the luminous SDSS DR7 quasars in the same redshift range, whose median FWHM is estimated to be 5266 km s −1 , the line width of the QOP=4 quasars is systematically smaller.Seen from the constant BH mass relation plotted by black dashed lines, the QOP=4 quasars are expected to harbor less-massive SMBHs than the SDSS quasars.We note that there is one quasar measured to have FWHM less than 1000 km s −1 .By visually inspecting the spectrum, we find that its C IV emission line can be fit with a single Gaussian model with FWHM< 1000 km s −1 .Thus, we classify this object as a narrow-line AGN, and reject it in the broad-line quasar sample.
The uncertainty of the continuum and emission line measurements is determined by fitting mock spectra generated from the best-fit model of each quasar.We add Gaussian random noise, which is determined by the standard deviation of the residuals of the fitting, to the best-fit model spectra.Then, we apply the same fitting procedure to each mock spectrum, and measure its line width and continuum luminosity.For each quasar, we construct 50 mock spectra.The 1σ uncertainty in the C IV line width and L 1350 measurements are then evaluated by the rms scatter of those measured values from the mock spectra.Table 7 summarizes the uncertainties of the 80 QOP=4 quasars.

BH mass, bolometric luminosity and Eddington ratio
Based on the FWHM of C IV 1549 Å emission line and continuum luminosity L 1350 of the 80 QOP=4 quasars at 3 < z < 4.5, we can estimate their virial BH mass following equation ( 6).The derived BH mass covers a range of 7.4 < log M BH /M ⊙ < 9.7, with a median of log M BH /M ⊙ = 8.5.For each quasar, the 1σ error in the virial BH mass is calculated from the rms scatter of the 50 generated mock spectra as described in section 3.2.Table 7 summarizes the results.
There are quasars whose C IV emission line shows an extended blue wing, which is fitted by a blueshifted broad component.It is known that such blueshift correlates with spectral properties in quasar spectrum, such as equivalent width (EW) of C IV emission line, suggesting the blueshifted component may associate with a nonvirial motion (e.g., Richards et al. 2011).Thus, the necessity to correct for the contribution of the non-virial component to the C IV emission line profile in the virial BH mass estimates has been claimed (e.g., Shen & Liu 2012).Coatman et al. (2017) find the virial BH mass can be over-estimated by five times when the C IV emission line is strongly blueshifted by 3000 km s −1 .Here, since the wavelength range covered by our spectrum does not include other strong emission lines to directly determine the systemic redshift, we do not consider the blueshift correction in the virial BH mass estimates of the z ∼ 4 quasars in baseline.How the blueshift affects the virial BH mass estimates will be further discussed in section 5.5.
We estimate the bolometric luminosity, L bol , from the monochromatic continuum luminosity at 1350 Å with L bol = κλL 1350 , where κ = 3.81 is the bolometric cor-rection factor determined by the composite quasar SED (Richards et al. 2006b).Here, we note that the bolometric correction may depend on luminosity (e.g., Netzer et al. 2007b;Lusso et al. 2012;Trakhtenbrot & Netzer 2012;Duras et al. 2020;Shen 2020).But it becomes weak for the optical bolometric correction, especially in the low luminosity ranges (e.g., Lusso et al. 2012;Runnoe et al. 2012;Shen 2020).Considering the relatively low luminosity range covered by our quasars, we neglect the luminosity dependence of the bolometric correction.The resulting bolometric luminosity of the 80 QOP=4 quasars spans a range of 45.41 < log L bol /erg s −1 < 46.93.
The Eddington ratio is then calculated with λ Edd = L bol /L Edd , where L Edd is the Eddington luminosity given by L Edd = 1.26 × 10 38 (M BH /M ⊙ ) erg s −1 .For the 80 quasars, the estimated Eddington ratios distribute in a range of −1.4 < log λ Edd < 0.5, with a median of log λ Edd = −0.4.For each quasar, the 1σ error in the Eddington ratio is calculated from the rms scatter of the 50 generated mock spectra described in section 3.2.Table 7 summarizes the results.

Comparison with z ∼ 4 quasars in literature
We compare the distribution of the estimated virial BH mass of the 80 QOP=4 quasars at 3 < z < 4.5 with literature.Shen et al. (2011) compile similar measurements on the C IV 1549 Å emission line width and its continuum luminosity for quasars at z > 2 in the SDSS DR7 quasar catalog.The virial BH masses are estimated using the same calibration with the C IV FWHM (i.e., equation 6).We use those estimates in the comparison.
In Figure 8, we compare the distributions of virial BH mass, Eddington ratio and bolometric luminosity of the QOP=4 quasars and SDSS DR7 quasars at 3 < z < 4.5.The QOP=4 quasars identified in this work covers the luminosity range down to log L bol /erg s −1 = 45.4,which is more than one order of magnitude lessluminous than the SDSS quasars in the same redshift range.Meanwhile, the QOP=4 quasars distribute in the BH mass range down to 10 7.4 M ⊙ with a median of log M BH /M ⊙ = 8.5, which is around one order of magnitude less-massive than the SDSS quasars, whose median BH mass is estimated to be log M BH /M ⊙ = 9.3.On the other hand, both of the quasar samples have the Eddington ratio similarly distributing around log λ Edd ∼ −0.5, with the QOP=4 quasars slightly extend towards the high Eddington ratios.The lower luminosities of the QOP=4 quasars are thus mainly driven by their smaller BH masses.In addition, we see a wider dispersion in virial BH mass for the QOP=4 quasars, which is consistent with previous studies focusing on deep AGN sam- ples at z < 2 (e.g., Merloni et al. 2010;Schulze et al. 2015).
In Figure 9, we show the observed bivariate distributions of virial BH mass, Eddington ratio and luminosity of the QOP=4 quasars and SDSS DR7 quasars at 3 < z < 4.5.As plotted in the left panel, there is a positive correlation between the luminosity and BH mass.The ratio between the two quantities corresponds to the Eddington ratio, and the distribution suggests that these quasars are active SMBHs accreting at λ Edd > 0.01.Meanwhile, we find there is a sharp decrease towards smaller Eddington ratios of λ Edd < 0.1.It can be caused by the flux limit, especially in the lessmassive range.As shown by the lower dotted line in the right panel, even the flux limit of less-luminous quasars is still too shallow to detect those population in the lessmassive range of < 10 9 M ⊙ .
While there are studies suggesting the C IV emission line profile can be used as a single-epoch virial BH mass estimator (e.g., Vestergaard & Peterson 2006), some others indicate a large scatter between the line width of the C IV and Hβ emission lines (e.g., Netzer et al. 2007a;Shen & Liu 2012), and claim the Hβ or Mg II emission line width and its corresponding continuum luminosity can be more reliable estimators.Furthermore, the C IV emission line can be affected by the non-virial motion as well as the narrow or broad absorption lines (for review see Shen 2013).Therefore, we further compare the virial BH mass of the 80 QOP=4 quasars to that of quasars estimated using the Hβ or Mg II emission line at z ∼ 3.5, which are collected from Shemmer et al. (2004), Netzer et al. (2007a), Zuo et al. (2015), Saito et al. (2016), Trakhtenbrot et al. (2016) and Schulze et al. (2018).Results are plotted in Figure 10.
The quasars identified in this work fall in a similar luminosity range as those in Netzer et al. (2007a) and Trakhtenbrot et al. (2016).Both of the latter two studies adopt the Hβ emission line as the virial BH mass estimator.We only pick up quasars at 3 < z < 4.25 from their samples for the comparison.Overall, the QOP=4 quasars distribute in the similar ranges of virial BH mass and Eddington radio to their samples, while our quasars extend towards the low-mass and high-Eddington-ratio regime.The higher Eddington ratio of the QOP=4 quasars could be partly due to the systematic uncertainty of BH mass estimates with the C IV emission line.For example, the C IV blueshift can result in an underestimated BH mass and thus over-estimated Eddington ratio for the less-luminous quasars (for detailed discussions see section 5.5).Regardless of the systematic effects, for the z ∼ 3.5 quasars whose BH masses are estimated with the reliable Hβ or Mg II emission line, they also show a sudden decrease towards log λ Edd < −1, similar to the QOP=4 quasars.In order to determine the z = 4 BHMF and ERDF, we limit our sample to the 52 that meet the c1 criteria (the primary selection window of the z ∼ 4 quasars in section 2.1), since the selection function of the c1 criteria is well evaluated as a function of magnitude and redshift (Akiyama et al. 2018).We refer the sample as the less- In addition, we construct a luminous quasar sample in the same redshift range from the SDSS DR7 quasar catalog.In the SDSS legacy survey, Richards et al. (2002) select quasar candidates by stellar morphology and multi-color criteria, and apply a uniform target selection for spectroscopy.Its selection function is evaluated as a function of redshift and magnitude by Richards et al. (2006a).At z ∼ 4, the selection efficiency is close to 100% down to i SDSS = 20.2.The SDSS DR7 catalog contains the full quasar sample (Schneider et al. 2010).As described in section 3.4, the virial BH masses of quasars in the catalog are estimated using the same C IV calibration in this work (i.e., equation 6) at z > 2 (Shen et al. 2011).
For the luminous quasar sample, firstly, we only consider quasars with "uniform flag" of 1, which means the quasars are uniformly selected by the final quasar selection algorithm in Richards et al. (2002).Then, we limit the redshift range of quasars to be 3.5 < z < 4.25 to match that of the less-luminous quasar sample.There are 2,144 quasars meeting the criteria.Finally, to prevent overlaps in the selection function on the UV ab-He et al.
solute magnitude M 1450 vs. redshift plane with that of the less-luminous quasar sample, we only select quasars above the magnitude limit of M 1450 = −25.75,which is determined by the upper limit of the selection function of the less-luminous quasar sample (for details see section 4.2).Here, we evaluate M 1450 of SDSS quasars by their virial BH masses and Eddington ratios provided by Shen et al. (2011).We firstly derive the bolometric luminosity through L bol = λ Edd L Edd .Then, we obtain M 1450 of the SDSS quasars by adopting the bolometric correction in Runnoe et al. (2012) as follows: log L bol = 4.74 + 0.91 log λL 1450 . (7) In total, we select 1,462 quasars from the SDSS DR7 quasar catalog as the luminous quasar sample.
The luminous quasar sample covers the BH mass range of 8.1 < log M BH /M ⊙ < 10.67, with a median of log M BH /M ⊙ = 9.49; and their Eddington ratios distribute in the range of −1.71 < log λ Edd < 0.81, with a median of log λ Edd = −0.49.Distributions of the BH mass, Eddington ratio and luminosity of the luminous quasar sample can be seen by grey filled histograms in Figure 8 and grey dots in Figure 9.We confirm the distributions of BH mass and Eddington ratio of the luminous sample are consistent with those of the entire SDSS DR7 quasars at 3.5 < z < 4.25.

Broad-line AGN BHMF and ERDF with the V max method
We firstly determine the broad-line AGN BHMF and ERDF at 3.5 < z < 4.25 using the V max method (the binned BHMFs and ERDFs hereafter; Avni&Bahcall 1980), which is applied to determine the AGN luminosity function in previous studies.The mass range of 7.5 ≤ log M BH /M ⊙ < 11 is divided into 12 bins with a bin width of ∆logM BH = 0.3 dex.The number density in the bin between (log where i is the index of each quasar, and n is the total number of quasars in the mass bin.Here, V i is the comoving effective survey volume of the i-th quasar in the mass bin, and 1/V i represents the contribution of the i-th quasar to the number density of the mass bin.
The effective survey volume between z min and z max is calculated with where d A is the angular diameter distance, c is the speed of light, and dτ /dz = 1/ [H 0 (1 + z)E(z)] is the lookback time at z. Here, E 2 (z) = Ω m × (1 + z) 3 + Ω Λ .Ω(M 1450,i , z) represents the effective survey area of the quasar with M 1450,i and z.The 1σ uncertainty of the number density is calculated following the Poisson statistics by In order to evaluate the effective survey area Ω(M 1450,i , z) of the less-luminous quasar sample, we adopt the effective survey area of the c1 quasars determined as a function of absolute magnitude and redshift in Akiyama et al. (2018) (see Fig. 16 therein).The less-luminous quasars are not uniformly selected for spectroscopy on the absolute magnitude vs. redshift plane.We thus estimate their survey area by multiplying the fraction of the spectroscopically-identified quasars among the entire c1 quasars in each absolute magnitude and redshift bin.We divide the absolute magnitude range of −27 ≤ M 1450 ≤ −20 into 70 bins with a bin width of ∆M 1450 = 0.1, and the redshift range of 3.05 ≤ z ≤ 4.85 into 18 bins with a bin width of ∆z = 0.1.In each bin (M 1450,i , z j ), we evaluate the fraction considering a range of (M 1450,i − 0.19) ∼ (M 1450,i + 0.19) and (z j − 0.19) ∼ (z j + 0.19).The range is chosen to ensure a sufficient number of c1 quasars in each bin while avoiding excessive smoothing.We use spectroscopic redshifts for spectroscopically-identified quasars, and photometric redshifts for the remaining targets.We do not correct for the fraction of the identified contaminants in each bin since it is negligibly small as described in section 2.3.The resulting effective survey area of the less-luminous quasar sample is shown with orange contours in the right panel of Figure 11.For the luminous quasar sample, we directly adopt the effective survey area of 6,248 deg 2 estimated by Shen & Kelly (2012).
The comoving effective survey volume for each quasar is obtained by integrating equation 9 over the redshift range between z min = 3.5 and z max = 4.25.While it is suggested the observed number density of AGNs rapidly evolves at z < 3 (e.g., Hasinger et al. 2005), the evolution becomes milder at z > 3 (e.g., Ueda et al. 2014).Since our quasar samples cover a narrow redshift range (3.5 < z < 4.25), we do not consider the redshift dependence in the integration for simplicity.We confirm the change would be negligibly small even if we introduce the redshift dependence.
The resulting binned broad-line AGN BHMFs determined by the less-luminous, luminous and combined quasar samples are plotted by orange open squares, grey open triangles and red circles in the upper left, upper middle and upper right panels of Figure 12, respectively.Tabel 8 summarizes the results.The BHMF of the luminous quasar sample peaks at log M BH /M ⊙ ∼ 9.6, but the turnover disappears in the BHMF of the combined quasar sample.Thus, the peak is likely due to incompleteness of the luminous quasar sample in the less-massive bins.The BHMF of the combined sample covers a wide mass range of 7.5 ≤ log M BH /M ⊙ ≤ 10.5 with the massive end dominated by the luminous sample and the less-massive end dominated by the less-luminous sample.There is a sharp decline in the high mass end above log M BH /M ⊙ > 9.6 as well as a mild decline in the low mass end below log M BH /M ⊙ < 8.5.The decline in the low mass end suggests the BHMF in the less-massive bins may be still affected by incompleteness even after combining the less-luminous quasar sample.
The binned broad-line AGN ERDFs are evaluated in the same way as for the binned broad-line AGN BHMFs.We divide the Eddington ratio range of −2 ≤ logλ Edd < 1 into 11 bins with an interval of ∆ log λ Edd = 0.3 dex.The number density in each bin is determined over the range of (log λ Edd − ∆ log λ Edd /2) ∼ (log λ Edd + ∆ log λ Edd /2).The 1σ uncertainty of the number density is calculated following the Poisson statistics.
The resulting binned broad-line AGN ERDFs determined by the less-luminous, luminous and combined quasar samples are plotted by orange open squares, grey open triangles and red circles in the lower left, lower middle and lower right panels of Figure 12, respectively.The results are summarized in Tabel 9.The ERDFs of the luminous and less-luminous quasar sam-ple cover the similar range of −1.5 ≤ log λ Edd ≤ 0.5.Both of them have similar shape with a turnover around log λ Edd = −0.5.Since the number density of lessluminous quasars is much higher than that of luminous quasars, the ERDF of the combined sample is dominated by the less-luminous quasars.

Broad-line AGN BHMF and ERDF with the Maximum Likelihood Method
Though we consider the effective survey area in determining the binned broad-line AGN BHMFs and ERDFs, the results can still be affected by the selection incompleteness caused by the flux limit of the quasar samples.As can be seen in Figure 9, with a certain flux limit, only quasars with sufficiently high Eddington ratios can be observable in the less-massive bins, and only quasars with large enough BH mass can be detected in the low-Eddington-ratio bins.Quasars with low BH masses and low Eddington-ratios can be missed as they are below the flux limit.The binned BHMFs and ERDFs constrained by using the flux-limited quasar samples can thus be incomplete in the less-massive and lower-Eddington-ratio ranges.
The BHMFs and ERDFs corrected for the flux limit (the intrinsic BHMFs and ERDFs hereafter) can be derived statistically by assuming multiple functional shapes for the intrinsic BHMF and ERDF (e.g., Schulze & Wisotzki 2010;Nobuta et al. 2012;Shen & Kelly 2012;Kelly & Shen 2013;Schulze et al. 2015).In this work, we follow the maximum likelihood method introduced in Schulze et al. (2015).The method enables to combine multiple samples with different selection functions, and He et al.
In Schulze et al. (2015), the intrinsic broad-line AGN BHMF and ERDF are determined jointly through the intrinsic bivariate distribution function Ψ(M BH , λ Edd , z), which can be treated as the multiplication of BHMF and ERDF component, i.e., Ψ(M BH , λ Edd , z) = ϕ * × ψ Edd × ψ BH , where ϕ * is the normalization.Here, a mass dependence in the ERDF component is allowed, i.e., ψ Edd = ψ Edd (M BH , λ Edd , z).The intrinsic ERDF ϕ Edd is then given by integrating the bivariate distribution function over log M BH between log M BH,min /M ⊙ = 7.5 and log M BH,max /M ⊙ = 11, i.e., The intrinsic BHMF ϕ BH is also given by the same integration over log λ Edd between log λ Edd,min = −2 and log λ Edd,max = 1.We assume the BHMF can be described either with a double power-law function or a Schechter function since the double power-law function can reproduce the observed AGN luminosity functions up to z ∼ 6, and the Schechter function can represent the observed stellar mass functions of galaxies.We assume the ERDF follows a log-normal function as the observed ERDFs at z > 3 in literature show clear turnover, similar to the log-normal distribution (e.g., Kollmeier et al. 2006;Willott et al. 2010;Shen & Kelly 2012;Shen 2019;Farina et al. 2022).Meanwhile, there are also studies indicating the turnover can be caused by shallow detection limits of quasar surveys, and the Schechter function may be more reliable in describing the intrinsic ERDF (e.g., Shen & Kelly 2012;Schulze et al. 2015;Li et al. 2023).Thus, we also adopt a Schechter function for the intrinsic ERDF The BH mass dependence of ERDF is taken into account by assuming the log λ ′ Edd term varies with BH mass in a linear form of log λ ′ Edd = log λ * Edd + k λ (log M BH − log M BH,l ), where log M BH,l = 7 is fixed.A second-order polynomial function is also considered in the mass dependence, but we find the best-fit secondorder coefficient is consistent with zero.We thus only consider the linear term.The redshift dependence of BHMF and ERDF is not included in the fitting since the quasar samples cover a narrow redshift range of 3.5 < z < 4.25.
There are two effects affecting the intrinsic distribution Ψ(M BH , λ Edd , z) to the observed one Ψ o (M BH , λ Edd , z): one is the uncertainties associated with the virial BH mass estimates and the bolometric correction, and another one is the selection function of quasar samples over the M BH -λ Edd plane.The former effect can be considered by convolving the intrinsic bivariate distribution function Ψ(M BH , λ Edd , z) with the uncertainties of virial BH mass estimates σ VM and of bolometric correction where the uncertainties σ VM and σ BC are assumed to follow a log-normal function as M BH,c and M BH are the virial and true BH mass, respectively.L bol,c and L bol are the converted and true bolometric luminosity, respectively.For the uncertainty in the virial BH mass estimates, we adopt σ VM = 0.2 dex following Schulze et al. (2015); for the uncertainty in the bolometric correction, we assume σ BC = 0.1 dex (Richards et al. 2006b;Shen et al. 2008).The latter effect can be evaluated by directly multiplying the convolved bivariate distribution function Ψ c (M BH,c , λ Edd,c , z) with the effective survey area Ω(M BH , λ Edd , z), i.e., Ψ o (M BH,c , λ Edd,c , z) = Ω(M BH,c , λ ,cEdd , z) The effective survey area of the luminous and lessluminous quasar samples is described as a function of absolute magnitude and redshift in section 4.2.To obtain the effective survey area as a function of BH mass and Eddington ratio, we convert Ω(M 1450 , z) to Ω(M BH , λ Edd , z).We follow the method in Kelly & Shen (2013).Firstly we calculate a bolometric luminosity L bol for any pair of BH mass M BH and Eddington ratio λ Edd , and then transform it to the monochronic luminosity at 1450 Å through the bolometric correction in Runnoe et al. (2012) (i.e., equation 7).We confirm the converted monochronic luminosity at 1450 Å is consistent with that directly measured from the spectrum.The effective survey area at the corresponding absolute magnitude M 1450 and redshift z is then adopted as the effective survey area at that BH mass, Eddington ratio and redshift.
For multiple quasar samples, if the respective survey area is independent with each other, their effective survey areas can be combined to one, i.e., Ω(M BH , λ Edd , z) = j Ω j (M BH , λ Edd , z), where j is the index of the respective quasar sample (Schulze et al. 2015).As mentioned in section 4.1, the luminous quasar sample is carefully selected not to overlap with the less-luminous quasar sample in absolute magnitudes.The respective survey area of the two quasar samples can thus be treated as independent area, and can be directly summed up for the combined quasar sample.In Figure 13, the comoving effective survey volume V eff of the two quasar samples is plotted against BH mass and Eddington ratio.With the observed bivariate distribution function of BH mass and Eddington ratio Ψ o (M BH , λ Edd , z), the probability of detecting the i-th observed quasar at its observed BH mass, Eddington ratio, and redshift can be written with where (20) is the normalization.We then apply the maximum likelihood method to minimize the likelihood function S = −2 ln L, where L = n i p i is the likelihood of detecting the observed quasar sample (Marshall et al. 1983), and n is the total number of quasars.The likelihood function can be written by The latter term can be simplified to 2n ln N if Ψ o is constant for all of the quasars.We minimize the likelihood function S by applying the downhill simplex algorithm (Nelder & Mead 1965).The free parameters are M * , α and β for the intrinsic BHMF; and λ * Edd , α λ /σ λ and k λ for the intrinsic ERDF.ϕ * is the normalization of the best-fit models, and it is derived by scaling the number count of quasars predicted by the best-fit models to that of the observed quasars, i.e., ϕ * = N obs /N model , where N obs =52, 1462 and 1514 for the less-luminous, luminous and combined quasar samples, respectively.N model can be obtained by integrating equation 20 over BH mass, Eddington ratio and redshift.The integration ranges are 7.5 < logM BH /M ⊙ < 11, −2 < logλ Edd < 1, and 3.5 < z < 4.3.For the luminous quasar sample, we further limit the integration range of BH mass to 8.5 < logM BH /M ⊙ < 11.
Uncertainty of the best-fit parameters is determined following the descriptions in Pawitan (2001).For the maximum likelihood estimates, the inverse of the observed Fisher information matrix can be the asymptotic covariance matrix, and the 1σ uncertainties of the bestfit parameters are then the square roots of the diagonal elements of the covariance matrix.Here, the observed Fisher information matrix is the second-order derivative of the log-likelihood evaluated at the maximum likelihood estimates, i.e., where θ ML are the best-fit parameters, and 1 ≤ i, j ≤ n para are the index of parameters.Uncertainty of the k-th parameter (1 ≤ k ≤ n para ) can then be calculated by The Fisher information matrix does not directly provide uncertainty of the entire best-fit models.Here, we determine that uncertainty following the χ 2 assumption.If the uncertainty of maximum likelihood estimates follows a χ 2 distribution, the 1σ uncertainty of the best-fit models can be determined when the likelihood function S increases from its best-fit value by 1.We then search for models having the change in S by one from its bestfit value.For each parameter, we randomly change it from its best-fit value and fix it.We then re-minimize the likelihood function S with the remaining parameters, and measure how much the minimum of S changes from its best-fit value.The same procedure is repeated for 50 times for each parameter, and all the models having the minimum of S within one from its best-fit value are adopted.The upper and lower boundaries of the adopted models are regarded as the 1σ uncertainty.We confirm the measured 1σ uncertainties are consistent with those evaluated by the Fisher information matrix.

The intrinsic broad-line AGN BHMFs and ERDFs
The intrinsic broad-line AGN BHMFs and ERDFs derived with the maximum likelihood method are shown in Figure 12.The functions determined with the lessluminous, luminous and combined quasar samples are plotted by orange, grey and red lines in the left, middle and right panels, respectively.Table 10 summarizes the best-fit parameters with the 1σ uncertainty of the respective model.
All of the four parametric models can get converged with the combined quasar sample, and they yield broadly consistent BHMFs.As plotted by lines in the left and middle panels of Figure 12, large discrepancies among the models only appear at the high-mass end for the less-luminous quasar sample, and at the lowmass end for the luminous quasar sample.In those BH mass ranges, the number of the observed quasars in the respective sample is limited, resulting in large scatter in different models.After combining the two quasar samples, that scatter disappears.For the ERDFs, there is still a difference among parametric models, especially between the log-normal and Schechter models, even after combining the quasar samples.The difference may be caused by the turnover functional shape of the lognormal function.The pink shaded areas in the figure represent the 1σ uncertainties of the double-power-law BHMF and log-normal ERDF model.These uncertainties only become significant at the low-mass and low-Eddington-ratio ends.
We compare the intrinsic BHMFs with the binned ones.The intrinsic BHMFs of the less-luminous and luminous quasar sample show broadly consistent number densities with the binned BHMFs over the mass ranges of 7.5 < log M BH /M ⊙ < 9.5 and 9.5 < log M BH /M ⊙ < 10.2, respectively.In the higher mass ranges, the intrinsic BHMFs decline more sharply than the binned ones; while in the lower mass ranges, the intrinsic BHMFs show higher number densities than the binned ones.The former is caused by correcting for the uncertainty associated with the virial BH mass estimates and the bolometric correction, and the latter is the consequence of accounting for the incompleteness of quasar samples due to their flux-limit selections.After combining the two quasar samples, the intrinsic and binned BHMFs show basically consistent number densities in the range of 8 < log M BH /M ⊙ < 10.2, and the flux-limit correction only becomes significant around log M BH /M ⊙ ∼ 7.5.For the ERDFs, in the combined quasar sample, the intrinsic log-normal ERDFs well follow the binned ones, while the intrinsic Schechter ERDFs suggest a flattened low-Eddington-ratio end at log λ Edd < −1 without clear turnover.
As described in section 4.3, in the maximum likelihood fitting, we consider the uncertainty in the virial BH mass estimates by assuming σ VM = 0.2 dex and that in the bolometric correction by adopting σ BC = 0.1 dex.It is suggested that the σ VM can reach ∼0.4 dex for an individual BH mass estimate with the C IV -based calibration (e.g., Vestergaard & Peterson 2006;Coatman et al. 2016).Here, we examine how these two uncertainties affect the determination of BHMFs by assuming multiple uncertainties of σ VM = 0, 0.1, 0.2, 0.3 and 0.4 dex, and σ BC = 0, 0.1, 0.3 dex.The resulting BHMFs and ERDFs are plotted by lines in Figure 14.
Compared to the binned BHMF, which does not consider the uncertainty, assuming a σ VM uncertainty results in a steepened massive end of the intrinsic BHMF.The high-Eddington-ratio end of intrinsic ERDF also becomes sharper, while the change is relatively small.On the other hand, assuming a σ BC uncertainty leads to narrow intrinsic ERDF.The less-massive end of intrinsic BHMF also becomes slightly steepened.Correcting for larger uncertainties can yield a more steepened highmass (high-Eddington-ratio) end in the intrinsic BHMF (ERDF), while the less-massive end does not show significant change.In this study, since we focus on how much the uncertainties can affect the distributions of the statistical sample instead of an individual estimate, we fix σ VM = 0.2 dex and σ BC = 0.1 dex hereafter, consistent with previous work (e.g., Kelly & Shen 2013;Schulze et al. 2015).
Thanks to the wide BH mass ranges covered by the combined quasar sample, we are able to investigate the mass dependence of ERDFs at z = 4.In Figure 15, we plot the intrinsic bivariate distribution function of BH mass and Eddington ratio.No matter the Schechter or log-normal model is applied for ERDF, the mass dependence of ERDF is visible.Massive quasars are likely to have lower Eddington ratios, while the less-massive ones are likely to have higher Eddington ratios.
In the parametric model, the mass dependence of ERDFs is described with the parameter k λ .As can be seen from the best-fit models in Table 10, no matter which parametric model combination is adopted, k λ keeps negative even with the uncertainty, suggesting a mass dependence of ERDFs is required.For comparison, we applied the maximum likelihood method with k λ fixed to zero, and the resulting double-power-law BHMF and log-normal ERDF are plotted by green solid lines in the right panels of Figure 12.Much larger corrections for the incompleteness at the less-massive and low-Eddington-ratio ends are required in this case.Considering those large corrections, the number density in the low-mass and low-Eddington-ratio ranges of log M BH /M ⊙ < 8.5 and log λ Edd < −1, which corresponds to the absolute magnitude range of M 1450 > −22, is predicted to exceed the observation constraints of ∼ 10 −6 Mpc −3 (e.g., Akiyama et al. 2018).We confirm that the predicted number density will not be largely changed if the Schechter function is applied to ERDF.
In order to further examine the mass dependence of ERDFs, we divide the combined quasar sample into the high mass bin of log M BH /M ⊙ > 9.5 and the low mass bin of log M BH /M ⊙ < 9.5, where log M BH /M ⊙ = 9.5 is the median BH mass of the combined quasar sample.In the respective mass bins, the binned ERDFs are evaluated by the same V max method described in section 4.2.Additionally, we acquire the intrinsic ERDFs in the two mass bins by integrating equation 11 over the corresponding mass ranges.Since the Schechter and double-power-law BHMFs are roughly consistent, we focus on the latter model hereafter for simplicity.Results are plotted in Figure 16.Except for the low-Eddington-ratio ends largely affected by the flux-limit corrections, the intrinsic mass-dependent ERDFs well follow the binned ones in both of the mass bins.By comparing the ERDFs in the two mass bins, the mass dependence of ERDFs can be clearly seen: the ERDFs in the high mass bin tend to peak at lower Eddington ratio of log λ Edd ∼ −1, while those in the low mass bin peak at higher Eddington ratio of log λ Edd ∼ −0.5.The semi-analytic model simulating the mass assembly of BHs at z > 4 in Piana et al. (2021) shows a similar trend that more massive BHs tend to accrete at lower Eddington ratios at z ∼ 4 − 5, as they switch to the gas-limited accretion phase at higher redshifts.

Goodness of the Maximum Likelihood fitting
Maximum likelihood fitting does not directly provide the goodness of fit.Thus, we evaluate the goodness of fit by comparing the best-fit bivariate distribution functions above the detection limits with the distributions of the two quasar samples on the M BH and λ Edd plane.As plotted in Figure 17, for the luminous quasar sample, we see good agreements between the observed and model distributions with both of the ERDF models; and for the less-luminous quasar sample, the observed and model distributions are still consistent with each other, but the model distribution with the log-normal ERDF follows better with the observed distributions than that with the Schechter ERDF.We quantitatively evaluate the goodness of fit by the two-dimensional Kolmogorov-Smirnov (2DKS) test on the M BH and λ Edd plane (Fasano & Franceschini 1987).The resulting probabilities exceed 5% for both of the models, i.e., we can not reject the statement that the observed and model distributions are drawn from the same parent distribution.
The quasar luminosity functions, which can be obtained by convolving the intrinsic BHMFs and ERDFs, is another test of the goodness of fit.For the two quasar samples, using the same V max method mentioned in section 4.2, we evaluate their binned luminosity functions in the range of −30 < M 1450 < −20, with a bin width of 0.5.We then convolve the best-fit intrinsic BHMFs and ERDFs of the combined quasar sample to derive the intrinsic luminosity functions of quasars over the same luminosity range.
The resulting binned and intrinsic luminosity functions of quasars are plotted in the left panel of Figure 18.The observed and model luminosity functions show good He et al.Furthermore, following Bongiorno et al. (2016), we use the Akaike information criterion (AIC; Akaike 1974) to compare the fitting quality of the respective parametric model.It is defined as AIC= S + 2K + 2K(N/(N − K − 1)), where S is the minimum likelihood of fitting, K is the number of free parameters in model, and N is the size of sample.We list the resulting relative AIC for each parametric model in Table 10.
For fitting the combined quasar sample, BHMFs in the double power-law model get better quality than those in the Schechter model, and ERDFs in the lognormal model achieve smaller AIC value than those in the Schechter model.The parametric model with BHMF in the double power-law function and ERDF in the log-normal function has the minimum AIC value, i.e., it can be regarded as the best-fit model.We note that the relative AICs can be enhanced by ≳ 50 if no mass dependence of ERDF is assumed in the fitting.
In summary, for discussions hereafter, we focus on the best-fit double-power-law BHMF with the log-normal ERDF, determined by the combined quasar sample with the mass-dependent ERDF.We also put the best-fit double-power-law BHMF with the Schechter ERDF as a comparison.

Comparison with previous work
We compare the derived broad-line AGN BHMFs and ERDFs with those in literature.Shen & Kelly (2012) and Kelly & Shen (2013) examine the broad-line AGN BHMFs and ERDFs in the redshift range of 0.3 < z < 5.They utilize a uniformly-selected sample of ∼58,000 SDSS DR7 quasars.At z > 1.9, the virial BH masses of the quasars are estimated with the same C IV -calibrated single-epoch estimator used in this work.In both of the studies, a Bayesian method with flexible models describing the underlying BHMFs and ERDFs is adopted to correct for the selection biases introduced by the flux limits, and luminosity-dependent scatter in the virial BH mass estimate is taken into account.The resulting BHMFs and ERDFs at z = 3.75 in Shen & Kelly (2012) and Kelly & Shen (2013) are shown by sky-blue and olive lines in Figure 19, respectively, together with those derived in this work.
In the left panel of the figure, we find a broad agreement between BHMFs in this work and theirs in the high-mass range of log M BH /M ⊙ ≳ 9.The luminous quasars mainly determine the constraints in that mass range, and the luminous quasar sample used in this work is constructed from the same SDSS DR7 quasar catalog.Although different methods, i.e., the Bayesian and Maximum likelihood methods, and the functional models are adopted, the consistent intrinsic BHMFs can be obtained.Discrepancy in BHMFs appears in the lowmass range of log M BH /M ⊙ ≲ 8.5, where the completeness of the SDSS quasar sample drops below 10% (Kelly & Shen 2013).The discrepancy can thus be caused by the uncertainty introduced by the large incompleteness correction of the SDSS quasar sample.In this work, by combining the luminous and less-luminous quasar samples, we can constrain the low-mass end of the z = 4 BHMF down to log M BH /M ⊙ ∼ 7.5 with the completeness more than 10%.
In the right panel of Figure 19, the ERDFs derived in this work are broadly consistent with those in Shen & Kelly (2012) and Kelly & Shen (2013), especially at the high-Eddington-ratio end.The ERDFs show discrepancy at log λ Edd < −1, which can also be caused by the incompleteness of the SDSS quasar sample in this range as indicated by the black dashed line.

Cosmological evolution of the broad-line AGN BHMFs and ERDFs
In order to clarify what drives the "down-sizing" trend seen in the AGN luminosity functions, we compare the z ∼ 4 broad-line AGN BHMFs and ERDFs with those at 0 < z < 6 in literature: Shen & Kelly (2012) and Kelly & Shen (2013) covering 0.3 < z < 5 as described above; Schulze & Wisotzki (2010) using 329 z < 0.3 quasars and Seyfert-1 galaxies drawn from the Hamburg/ESO Survey; Schulze et al. (2015) with the 1.1 < z < 2.1 quasars from the VVDS epoch-2 AGN, the zCOSMOS 20k AGN, and the SDSS DR7 QSO samples; Nobuta et al. (2012) with the X-ray-selected 1.18 < z < 1.68 quasars in the Subaru XMM-Newton Deep Survey (SXDS); Willott et al. (2010)  z ≳ 6 quasars from the Canada-France High-z Quasar Survey (CFHQS); and Wu et al. (2022) with 47 quasars at 5.7 ≤ z ≤ 6.5 from SDSS.We note that different single-epoch virial BH mass estimators are adopted in each work: Shen & Kelly (2012) and Kelly & Shen (2013) use Hβ, Mg II and C IV emission lines at z < 0.7, 0.7 < z < 1.9 and z > 1.9, respectively; Schulze & Wisotzki (2010) use Hβ emission line; and other studies mostly rely on Mg II emission line.In addition, the z = 6 BHMF in Willott et al. (2010) already corrects for the obscured fraction of AGNs, including the Compton-thick ones, and the active fraction.It thus should be treated as an upper limit in the comparison with other broad-line AGN BHMFs, especially in the low-mass range as the obscuration correction becomes significant (e.g., Merloni et al. 2014).
In Figure 20, the cosmic evolution of the number density of the broad-line AGNs at a constant BH mass and Eddington ratio is shown.To guide the comparison, we connect the results at z ∼ 0 in Schulze & Wisotzki (2010), at z ∼ 1 − 2 in Schulze et al. (2015) and at z ∼ 4 in this work by dotted lines.Results from other studies broadly follow the lines.There is a large scatter in the lowest BH mass and Eddington ratio bins, which can be caused by large incompleteness corrections of the SDSS quasars in the ranges (Shen & Kelly 2012;Kelly & Shen 2013).In addition, we see the number density at constant Eddington ratios constrained by Kelly & Shen (2013) suddenly increases at z ∼ 2. Since the virial BH mass estimator is changed from Mg II to C IV at z ∼ 2 in their work, that sharp increase can be a consequence of the inconsistency between different virial BH mass estimators.Furthermore, results in Nobuta et al. (2012) show large excesses in the high Eddington ratio bins.It should be noted that the quasar sample in their work covers a relatively less-massive range of log M BH < 10, and they do not correct for the scatter induced by the virial BH mass estimate.The resulting BHMF and ERDF in their work thus show flattened massive and high-Eddington-ratio ends.
Following the dotted lines in the figure, from z = 2 to z = 6, the number density of massive SMBHs with log M BH /M ⊙ = 9.8 keeps constant up to z = 4 and then drops towards z = 6, while that of less-massive SMBHs with log M BH /M ⊙ ≲ 9 shows continuous declines towards z = 6.The number density of rapidlyaccreting SMBHs with log λ Edd ≳ −0.1 significantly increases from z ∼ 2 to z ∼ 4, while that of SMBHs with lower Eddington ratios of log λ Edd ≲ −0.5 decreases from z ∼ 2 to 4. Both of the evolutionary trends suggest the abundance of more massive and rapidly-accreting SMBHs peaks at higher redshifts, consistent with the "down-sizing" scenario.There are also semi-analytical simulations suggesting the number density of rapidlyaccreting SMBHs with log λ Edd ∼ 0 tends to peak at high redshifts (e.g., Shirakata et al. 2019).The possible cause of the trend may be due to the gas-rich environment of galaxies at high redshifts, making rapid accretion to easily occur.
In the less-massive range of log M BH /M ⊙ ≲ 9, we find the number density basically evolves with a similar slope from z = 2 to 4. The "up-sizing" evolutionary trend, which indicates a milder decline in the number density of less-massive or lower-Eddington-ratio SMBHs after the peaks towards high redshifts, is not seen in the comparison.

Obscuration correction and the active fraction of SMBHs at z ∼ 4
The current z ∼ 4 quasar sample only covers unobscured AGNs, here we estimate the BHMFs and ERDFs for the total AGN population by correcting for the obscured AGNs.Studies on hard-X-ray (2-10 keV) selected AGNs suggest the fraction of obscured AGNs among the entire AGN population deceases with increasing luminosity (e.g., Ueda et al. 2003;Hasinger et al. 2008;Merloni et al. 2014;Ueda et al. 2014).Physical origin of the dependence is not disclosed yet, and a possible scenario is that the inner edge of the dusty obscuring torus recedes with increasing luminosity (e.g., Lawrence 1991;Simpson 2005;Ricci et al. 2013;Toba et al. 2014;Vijarnwannaluk et al. 2022).Merloni et al. (2014) determine the relationship between the X-ray luminosity at 2-10 keV and the fraction of obscured AGNs among the entire AGN population with a complete sample of 1310 X-ray-selected AGNs at 0.3 < z < 3.5 in the XMM-COSMOS field.Here, we assume the relationship does not evolve with redshift, and adopt it as follows: Additionally, since the relationship is poorly constrained in the less-luminous range, we apply an upper limit of f obs as f max = 0.985, which is the obscured fraction at the luminosity of L 2−10 keV = 10 42 erg s −1 .It should be noted that the contribution of Compton-thick AGNs is not considered here since it is still highly uncertain.
The correction for the obscured fraction is applied in each BH mass and Eddington ratio bin on the M BH vs. λ Edd plane.Firstly we calculate the bolometric luminosity L bol of the combination of BH mass M BH and Eddington ratio λ Edd , and then convert the luminosity to the X-ray luminosity in the 2-10 keV band using the bolometric correction of L 2−10 keV = L bol /κ x , where the obscuration-corrected bivariate distribution against the Eddington ratio over −2 < log λ Edd < 1 and the BH mass over 7.5 < log M BH /M ⊙ < 11, respectively.
The resulting total active BHMFs and ERDFs are plotted with black lines in Figure 21.Solid and dashed lines show the results derived with the log-normal and Schechter ERDF, respectively.As expressed in equation 24, the obscured AGNs are assumed to dominate the less-luminous AGN population.We thus see significant correction for the obscured AGNs in the lowmass or low-Eddington-ratio ranges of the total active BHMFs and ERDFs.At the high-mass or high-Eddington-ratio ends, the contribution from the obscured AGNs can be negligibly small.We note that the correction in the model with lognormal ERDF is weaker than that with Schechter ERDF.The reason can be due to the sharp decline of the log-normal ERDF towards the low-Eddingtonratio end, where the obscured AGNs can dominate.Therefore, even after correcting for the obscured AGNs, the total spatial number density of active AGNs in the low-mass and low-Eddington-ratio regime predicted by the log-normal ERDF remains small.
Utilizing the total active BHMFs, we evaluate the active fraction of SMBHs at z = 4, which is defined as the fraction of active SMBHs with λ Edd > 0.01 among the total SMBH population, by comparing the total active BHMFs to the total BHMFs including quiescent SMBHs with λ Edd < 0.01.In order to derive the total BHMF, following Schulze et al. (2015), we convolve the total stellar mass function of galaxies Ψ * with the M BH − M star relation, i.e., where s = log M star − 11, and α, β and σ are the the normalization, slope and intrinsic scatter of the M BH − M star relation.
To compare with Schulze et al. (2015), we adopt the same M BH −M Bulge relation in McConnell & Ma (2013), i.e., α = 8.46, β = 1.05 and σ = 0.34, and the total stellar mass function of galaxies at z ∼ 4 in Ilbert et al. (2013).Here, as explained in Schulze et al. (2015), due to multiple uncertainties, such as the ratio between the stellar and spheroid mass at high redshifts, we use total stellar mass as a proxy of spheroid mass.
The M BH −M Bulge relation given by McConnell & Ma (2013) is examined in the local universe.Whether this relation evolves with redshift still remains controversial.In the high-redshift galaxies, the decomposition of their bulge components is not applicable.The stellar mass is then treated as a surrogate of bulge mass, and the M BH − M star relation is used instead (e.g., Reines & Volonteri 2015;Shankar et al. 2016;Suh et al. 2020).There are studies suggesting a possible redshift evolution in the M BH − M star relation (e.g., McLure et al. 2006;Decarli et al. 2010), while others suggest no evolution (e.g., Willott et al. 2017;Schulze & Wisotzki 2014;Suh et al. 2020).
Considering the uncertainties, in addition to the M BH −M Bulge relation given by McConnell & Ma (2013), we apply the M BH − M star relations determined in Suh et al. (2020) and Shankar et al. (2016).Suh et al. (2020) determine the relation for a sample of 100 X-ray-selected moderate-luminosity broad-line AGNs out to z ∼ 2.5 in the Chandra-COSMOS Legacy Survey, with the local AGNs in Reines & Volonteri (2015).The best-fit relation shows a steeper decline towards the low-M star end than that in McConnell & Ma (2013).Additionally, they also derive the relation with only the high-z AGNs.The resulting relation is consistent with that in McConnell & Ma (2013), but with a larger intrinsic scatter of ∼ 0.5 dex.Shankar et al. (2016) suggest the M BH − M star relation can be heavily biased by up to 50 times at small stellar masses due to the selection effect.We use the de-biased relation provided by Shankar et al. (2016).(2016) show more flattened low-mass ends with the normalizations over one order of magnitude smaller than the former two results.
The reason of this discrepancy can be caused by the sharp decline of the latter two relations in the low-mass range of log M star /M ⊙ < 11.Meanwhile, the massive end of the resulting total BHMF shows strong decline if the adopted M BH − M Bulge or M BH − M star relation has small intrinsic scatter.The results given with the relations in Shankar et al. (2016) and McConnell & Ma (2013) thus have a high-mass end declining even sharper than the broad-line AGN BHMFs.
With the total BHMFs, we evaluate the active fractions of broad-line AGNs as a function of BH mass.The results are presented by thin lines in the middle panel of Figure 22.The active fraction of broad-line AGNs shows a positive dependence on the BH mass: the lessmassive SMBHs with 7.5 < log M BH /M ⊙ < 8.5 show an average fraction of 0.01 ∼ 0.3, while the massive ones with log M BH /M ⊙ > 9 occupy 0.04 ∼ 0.8.Such BH mass dependence of the active fraction is also found in semi-analytical simulations of galaxy and BH evolution at z > 2 (e.g., Shirakata et al. 2020).
Beyond the peak around log M BH /M ⊙ = 9.5, the active fractions start decreasing.It should be noted that most of the applied M BH − M Bulge or M BH − M star relations are extrapolated in the mass range of log M BH /M ⊙ > 9.The converted total BHMF and the active fraction in that mass range thus become highly uncertain, and some resulting active fractions exceed unity.
The fraction of dark matter halos with a broad-line AGN is also examined by the clustering analysis at z ∼ 4. The estimated active fractions are 0.001 ∼ 0.06 for the less-luminous quasars (He et al. 2018), and 0.03 ∼ 0.6 for the luminous quasars (Shen et al. 2007).Although the methods to derive the fraction are different, we see broad agreements between the active fractions constrained in this study and those determined in the clustering analysis.The luminous SDSS quasars are expected to associate with massive SMBHs of log M BH /M ⊙ > 9 with high active fractions, while the less-luminous quasars are more likely to represent smaller SMBHs with lower active fractions.
With the correction for the obscured fraction, the active fractions of the total AGN population are plotted by thin lines in the right panel of Figure 22.They show constant values with a mean at 0.01 ∼ 0.15 over a wide BH mass range of 7.5 < log M BH /M ⊙ < 9, implying the phase of active accretion is occurring throughout the less-massive SMBH population without preferences on the BH mass.That active fractions are consistent with the values evaluated with the X-ray-selected AGNs in the deep survey fields (e.g., Bongiorno et al. 2012Bongiorno et al. , 2016;;Georgakakis et al. 2017), and with the semi-analytical model for the BH formation and growth (e.g., Li et al. 2023).Again, the increase beyond log M BH /M ⊙ > 9 is highly uncertain because the total BHMFs are evaluated with the extrapolated M BH − M Bulge or M BH − M star relations.
We investigate the cosmic evolution of the active fraction of the broad-line and total active AGNs by comparing the fractions at z = 4 with that at z = 2, which is estimated by Schulze et al. (2015) with the same method described above.Since only the M BH − M Bulge relation in McConnell & Ma (2013) is adopted in Schulze et al. (2015), we consider the result derived with the same relation at z = 4.While the uncertainty remains large, the massive SMBHs of log M BH /M ⊙ > 9.5 are likely to keep a high active fraction of ∼ 10% since the cosmic noon.We note that the comparison highly depends on the stellar mass function of galaxies and M BH − M Bulge relation at high redshifts.Future studies on examining the M BH − M Bulge relation over wide BH mass and redshift ranges can be important to improve confining the active fractions.We further examine the time evolution of BHMF between z = 4 and 6 using the intrinsic broad-line AGN BHMF and ERDF at z = 4 determined in section 4.4 with the total BHMF at z = 6 given by Willott et al. (2010).We try to reproduce the observed broad-line He et al.AGN luminosity functions at z = 4, 5, and 6, by optimizing the time evolution model of the broad-line AGN BHMF.We assume the evolution is driven only by the mass accretion during the AGN phase, and the evolution is described by the continuity equation (e.g., Yu & Tremaine 2002;Shankar et al. 2013): Here, ⟨ ṀBH ⟩ is the average accretion rate of SMBHs with BH mass M BH at time t.Assuming the SMBH growth is driven only by the AGN activity, this rate corresponds to the total mass accreted onto the SMBHs during the AGN phase, i.e., Ṁ BH = (1 − ϵ) Ṁin , where ϵ is the radiative efficiency and Ṁin is the total mass inflow rate.The remaining inflow mass is converted into luminosity L = ϵ Ṁin c 2 .As described in section 3.3, the Eddington ratio λ Edd is defined as the ratio between luminosity L and the Eddington luminosity L Edd = ℓM BH , where ℓ = 1.26 × 10 38 M −1 ⊙ erg s −1 .Thus, the average accretion rate can be written as where U and ⟨λ Edd ⟩ are the active fraction and the average Eddington ratio of SMBHs with BH mass M BH at time t, respectively.
By combining equation 28 and equation 29, the continuity equation can be modified to: Since solving the continuity equation backwards in time can easily result in unphysical BHMFs with negative number densities (e.g., Shankar et al. 2013), we start from z = 6, and then move forwards in time to z = 4 with a time step of ∆z = 0.1 to solve equation 30.
We use the total BHMF at z = 6 given by Willott et al. (2010) as the initial condition.At each time step, we adopt the total BHMF Φ(M BH , z) obtained from the previous step or from the assumed initial condition in the first time step.Then, the active fraction U is derived as the ratio between the active and total BHMF at the time step, i.e., U (M BH , z) = Φ AGN (M BH , z)/Φ(M BH , z).Here, we confirm the time step of ∆z = 0.1 is sufficiently small that smaller time steps yield consistent results.The total number of black holes between the two time steps is conserved.
At the time step, the active BHMF Φ AGN (M BH , z) can be estimated from the broad-line AGN BHMF ϕ BH (M BH , z) after correcting for the obscured fraction of AGNs, as described in section 5.3.We set the broad-line AGN BHMF ϕ BH to be free.Since the Schechter and double-power-law broad-line AGN BHMFs are consistent with each other at z = 4, we only consider ϕ BH with the Schechter function (equation 13) during z = 4 ∼ 6 to reduce the number of free parameters.We describe the time evolution of ϕ BH through the characteristic BH mass, i.e., log M * new = log M * +k MBH,1 (z−z 0 ); the slope, i.e., α new = α + k MBH,2 (z − z 0 ); and the normalization, i.e., log ϕ * new = log ϕ * + k MBH,3 (z − z 0 ) + k MBH,4 (z − z 0 ) 2 , where z 0 = 4.There are six free parameters in total, i.e., M * , α, k MBH,1 , k MBH,2 , k MBH,3 and k MBH,4 .ϕ * is the normalization at z = 4 evaluated in section 4.4.
Meanwhile, at the time step, we evaluate the average Eddington ratio in each BH mass and redshift through ⟨λ Based on the broad-line AGN ERDFs at z = 4 derived in section 4.4, we assume two models: the first one assumes no redshift evolution of ERDF during z = 4 ∼ 6; and the second one allows for a time evolution of ERDF by adding an additional term to the characteristic Eddington ratio, i.e., log λ , where z 0 = 4.We only adopt the log-normal ERDF as it shows better fitting quality than the Schechter ERDF in Section 4.5 (see Table 10), and the shape of the observed ERDFs at z ∼ 6 can be better described by the lognormal function (e.g., Willott et al. 2010;Shen 2019;Farina et al. 2022).The coefficient is set to k λ,z = 0.11 so that the broad-line AGN ERDF can peak at the Eddington limit log λ Edd ∼ 0 at z = 6, which is consistent with Willott et al. (2010).We confirm using the Schechter ERDF does not significantly change the resulting best-fit broad-line AGN BHMF.
Finally, with the radiative efficiency ϵ and the six free parameters of ϕ BH as inputs, we can estimate the BH mass growth during the time interval through ⟨ Ṁ BH ⟩, and derive the total BHMF at the time step by integrating both sides of equation 30.The broad-line AGN BHMF is also the output.It should be noted that we adopt the total BHMF to be the active BHMF in the mass ranges where the latter exceeds the former to ensure the active fraction to be less than one.The broadline AGN BHMF over those mass ranges will also be modified accordingly by considering the obscured fraction of AGNs.Furthermore, we do not consider the time evolution of obscured fraction of AGNs during the period as it still remains uncertain.
For simplicity, we fix the radiative efficiency to ϵ = 0.05, 0.1 and 0.15.The six free parameters of ϕ BH are then determined by fitting the convolved luminosity function of broad-line quasars, i.e., Φ QSO (L, z) = d log λ Edd P QSO (λ Edd | M BH , z)ϕ BH (M BH , z), to the observations at z = 4 from Akiyama et al. (2018), z = 5 from Niida et al. (2020) and z = 6 from Matsuoka et al. (2018), respectively.These observed luminosity functions are evaluated over a wide luminosity range so that the shape of ϕ BH can be well determined over a wide mass range.The fitting procedure utilizes the least-square method, and the range of fitting is set to −28.5 < M 1450 < −21.5, where the observed luminosity functions are constrained with good statistics.
With each of the radiative efficiency ϵ of 0.05, 0.1 and 0.15, we adopt both of the ERDF models.The resulting luminosity functions of broad-line quasars at z = 4, 5 and 6 are plotted by red, orange and blue lines in Figure 23, respectively.We see the time-evolving ERDF model can yield higher number densities than the constant ERDF model at the luminous ends.The discrepancy may be due to the higher average Eddington ratios allowed by the former model, and more luminous quasars are enabled to be produced.
We compare the model luminosity functions with the observational ones at z = 4, 5 and 6: the predicted number densities at the faint ends are consistent with observations at all redshifts and radiative efficiencies, while those at the luminous ends tend to be underestimated by over one order of magnitude compared to the observations at large radiative efficiencies ϵ ≳ 0.1.When adopting the radiative efficiency of ϵ = 0.15, the discrepancies between the model and observed luminosity functions at z = 4 and 5 can be larger than three times at the luminous bins, even considering the time-evolving ERDF model.Thus, such large radiative efficiencies are likely to be discarded in our models.
To fully reproduce the observed luminosity functions at z = 4, 5 and 6 over −28.5 < M 1450 < −21.5, either a small radiative efficiency of ϵ ≲ 0.05 with the constant ERDF model or a moderate radiative efficiency of ϵ ≲ 0.1 with the time-evolving ERDF model is required.Such preference of small radiative efficiencies of ∼ 0.05 in the SMBH growth at z ≳ 4 is also stated, and explained as the result of more common super-Eddington accretion at higher redshifts in semianalytical simulations of galaxy and BH evolution (e.g., Shirakata et al. 2020).We note that this small radiative efficiency is close to the lower boundary for a geometrically thin disk (Shakura & Sunyaev 1973).It may indicate "radiatively inefficient accretion flows (RIAFs)" caused by low-Eddington (log λ Edd < −2) or super-Eddington (log λ Edd > 0) accretion (for review see Inayoshi et al. 2020).
In Figure 24, we plot the resulting total and broadline AGN BHMFs in the top and bottom panels, respec- tively.The results at z = 4, 5 and 6 are plotted by red, orange and blue lines, respectively.We see the results show a dependence on the radiative efficiency.With large radiative efficiencies ϵ ≳ 0.1, there is little evolution in the total BHMF from z = 6.Thus, the massive ends of the broad-line AGN BHMFs at z = 4 and 5 need to sharply decline so that they would not exceed the total BHMF at the respective redshift.The resulting low number density of massive SMBHs of log M BH /M ⊙ > 9 causes the lack of luminous quasars, so the model luminosity functions can be highly underestimated at the luminous ends.
With smaller radiative efficiencies ϵ < 0.1, the total BHMF can be significantly evolved from z = 6, especially at the massive end of log M BH /M ⊙ > 9.5.Compared with the total BHMFs convolved from the stellar mass function at z = 4, which are obtained in section 5.3, the number densities at the massive end of log M BH /M ⊙ > 9.5 can be well reproduced, while those at the less-massive end of log M BH /M ⊙ < 9 require even smaller radiative efficiencies of ϵ ≲ 0.05.
Under different small radiative efficiencies ϵ < 0.1, the resulting broad-line AGN BHMFs keep consistent at all redshifts, and the result at z = 4 is in accordance with the intrinsic broad-line AGN BHMF derived in section 4.4.We note that the time-evolving ERDF model can result in slightly higher BHMFs than the constant ERDF model due to the larger average Eddington ratios allowed by the former model.The best-fit parameters with the radiative efficiency of 0.05 are summarized in Table 11.
As both the total and broad-line AGN BHMFs at z = 4 can be best reproduced with the radiative efficiency of 0.05, we further examine the time evolution of broad-line AGN BHMF during z = 4 ∼ 6 under the radiative efficiency of 0.05.As shown in the bottom right panels of Figure 24, overall, the predicted broadline AGN BHMFs at z = 4, 5 and 6 keep similar shapes.
Only the normalization shows a major evolution.The similar pure density evolution trend can be also seen in the luminosity functions of quasars across z = 4 − 6, which are displayed in Figure 23.We note that both the constant and time-evolving ERDF models can result in the same evolution trend.

Effects of the C IV -based BH mass estimator
As mentioned in section 4.4, the C IV -based virial BH mass estimates can be under-or over-estimated due to the blueshifted component of the C IV emission line.Here, we evaluate its effect in estimating the virial BH masses using the relation between the C IV blueshift and the ratio of virial BH mass estimated by C IV and Hα lines, which is given by Coatman et al. (2017).We set an upper limit of 4000 km s −1 for the blueshift since the quasars used to determine the relation in Coatman et al. (2017) barely cover the range beyond the limit.
The wavelength coverage of the quasar spectra obtained in this work does not cover other emission lines, such as Mg II , to determine the systematic redshift.We thus adopt the empirical relation between the C IV FWHM and its blueshift.The relation is determined by fitting a second-order polynomial to 19 luminous SDSS quasars at 2 < z < 2.7 with both of the quantities available in Coatman et al. (2016).To prevent overcorrections, we limit the blueshift to be positive.Then, we derive the corrected BH mass with the C IV FWHM.
After correcting for the effect caused by the C IV blueshifts, the luminous and less-luminous quasar samples have the median BH mass of log M BH /M ⊙ = 9.4 and 8.7, respectively.Compared to the estimates without the blueshift correction, both of the quasar samples cover narrower mass range, and almost all massive SMBHs with log M BH /M ⊙ ≳ 9.6 in the luminous quasar sample are corrected to have smaller mass around log M BH /M ⊙ ∼ 9.5.Meanwhile, the median Eddington ratio of the less-luminous quasar sample drops to log λ Edd = −0.7,while that of the luminous quasar sample increases to log λ Edd = −0.4.
Considering the large uncertainty of the method described above, we further examine the effect by applying the calibration of the C IV -based single-epoch virial BH mass estimate in Park et al. (2017).It is suggested that the calibration can produce a similar BH mass distribution to that computed from the blueshift-corrected formula in Coatman et al. (2017) The systematic uncertainty of the calibration is evaluated to be 0.37 dex (Park et al. 2017).
In Figure 25, we plot the luminous and less-luminous quasar samples with the virial BH mass estimated by the calibration in Park et al. (2017) by magenta dots and open stars, respectively.Compared to the results derived with the calibration in Vestergaard & Peterson (2006), we find that both of the quasar samples show much narrower mass ranges when adopting the calibration in Park et al. (2017): the less-luminous quasar sample spans the range of 8.0 < log M BH /M ⊙ < 8.9 with a median of log M BH /M ⊙ = 8.4, and the luminous quasar sample covers the range of 8.6 < log M BH /M ⊙ < 9.6 He et al. with a median of log M BH /M ⊙ = 9.0.Meanwhile, the Eddington ratio range of the less-luminous quasar sample decreases to −0.8 < log λ Edd < −0.2 with a median of log λ Edd = −0.5, while that of the luminous quasar sample increases to −0.5 < log λ Edd < 0.5 with a median of log λ Edd = −0.0.
With the virial BH mass estimated with the calibration in Park et al. (2017), we derive the binned and intrinsic BHMFs and ERDFs at z = 4 for the combined sample with the same method introduced in section 4.3.The resulting BHMFs and ERDFs are plotted by magenta lines in the top and bottom panels of Figure 26, respectively.The best-fit parameters are summarized in Table 10.
The broad-line AGN BHMFs derived with the calibration in Park et al. (2017) show broad consistency with those obtained with the calibration in Vestergaard & Peterson (2006)  by magenta open diamonds, when adopting the calibration in Park et al. (2017), the number density of massive SMBHs with log M BH /M ⊙ = 9.8 significantly drops from z = 2 to 4, while that of less-massive SMBHs with log M BH /M ⊙ ≲ 9 mildly declines.Additionally, when adopting the calibration in Park et al. (2017), the mass dependence of ERDFs becomes positive with the mass-dependent term k λ significantly greater than zero, suggesting more massive SMBHs tend to accrete at higher Eddington ratios.The trend is consistent with the results determined with the Mg II -based single-epoch virial BH mass estimate at lower redshifts z = 1 − 2 (e.g., Kelly & Shen 2013;Schulze et al. 2015).
The main discrepancy between the calibrations in Park et al. (2017) and Vestergaard & Peterson (2006) is caused by the relaxed constraints of slope parameters in Park et al. (2017).Park et al. (2013) claim that treating the slope parameters in calibration as free parameters can mitigate the FWHM-dependent biases, which may be partly induced by the blueshifted component of the C IV emission line.However, we should also note that the relaxation of slope parameters in Park et al. (2017) may contradict the virial theorem (i.e., M BH ∝ FWHM 2 ).In addition, the AGN sample utilized to calibrate the virial BH mass estimates in their work does not cover massive SMBHs with log M BH /M ⊙ > 9, which may also result in large uncertainties in the virial BH mass estimates of luminous quasars.

SUMMARY
In this paper, we estimate the C IV -based virial BH mass of less-luminous quasars at z ∼ 4. The quasar candidates are selected down to i = 23.2 with the gband dropout colors and stellar morphology from the HSC-SSP S16A-Wide2 dataset.In total, we identify 151 quasars at z ∼ 4, and estimate the virial BH mass of 80 less-luminous quasars at 3.04 ≤ z ≤ 4.44.Our main results can be summarized as follows: 1.The less-luminous quasars span a BH mass range of 7.4 < log M BH /M ⊙ < 9.7 with a median of log M BH /M ⊙ = 8.5, which is around one order of magnitude less-massive than the luminous SDSS DR7 quasars in the same redshift range.
2. Both the luminous and less-luminous quasars have the Eddington ratio similarly distributed around log λ Edd ∼ −0.5, but the latter ones slightly extend towards the high Eddington ratios.
Then, we determine the broad-line AGN BHMFs and ERDFs at z = 4 with 52 less-luminous quasars covering the BH mass range of 7.4 < log M BH /M ⊙ < 9.7 at 3.50 ≤ z ≤ 4.25.To improve the constraints over a wide BH mass range, we further construct a sample of 1,462 luminous quasars in the same redshift range from the SDSS DR7 quasar catalog.The combined sample can cover a BH mass range of 7.4 < log M BH /M ⊙ < 10.7.We adopt both the V max and maximum likelihood method.The former method yields the binned results that only quasars above the flux limit are considered, while the latter one can derive the intrinsic results by accounting for the incompleteness in the low-mass and/or low-Eddington-ratio ranges caused by the flux limit.Our main results can be summarized as follows: 1.For the intrinsic broad-line AGN BHMFs and ERDFs, the parametric model with BHMF in dou-ble power-law function and ERDF in log-normal function shows the best-fit result.In the case, the intrinsic and binned BHMFs and ERDFs are broadly consistent with each other.Large discrepancies between the two only appear in the lessmassive range of log M BH /M ⊙ ∼ 7.5, where the completeness of the combined quasar sample drops below 50%.
2. A negative mass dependence of ERDF, i.e., the massive SMBHs accrete with low Eddington ratios, is suggested by the best-fit model.Without that mass dependence, the flattened faint-end of the observed luminosity function of quasars at z = 4 can not be reproduced.
3. The intrinsic broad-line AGN BHMF and ERDF follow the "down-sizing" evolutionary trend in comparison with those at the cosmic noon z ∼ 2.
The number density of the most massive SMBHs of log M BH /M ⊙ ∼ 10 roughly keeps constant from z = 2 to 4, while that of smaller SMBHs continuously declines during the period.Additionally, the number density of rapidly-accreting SMBHs of log λ Edd ≳ 0 is significantly enhanced since z ∼ 2, while that of SMBHs with lower Eddington ratios continuously declines from z = 2 to 4. The abundance of massive and/or rapidly-accreting SMBHs tends to peak at high redshifts.
Furthermore, we correct for the fraction of the obscured AGNs to estimate the BHMF and ERDF for the total AGN population at z = 4.We then compare the intrinsic broad-line and total active AGN BHMFs to the total BHMF at z = 4.The results correspond to the active fractions of broad-line and total active AGNs among the entire SMBH population as a function of BH mass at z = 4.In order to obtain the total z = 4 BHMF, we convolve the total stellar mass function of galaxies with the M BH − M * relation.The main findings can be summarized as follows: 1.While uncertainties remain large, the active fraction of broad-line AGNs suggests a positive dependence on the BH mass: the less-massive SMBHs with 7.5 < log M BH /M ⊙ < 8.5 show an average fraction of 0.01 ∼ 0.3, while the massive ones with log M BH /M ⊙ > 9 occupy 0.04 ∼ 0.8. 2. The best-fit time evolution model suggests the broad-line AGN BHMF basically evolves in parallel during z = 4 ∼ 6 without significant changes in the shape.Meanwhile, the average Eddington ratios at z = 5 and 6 tend to increase compared to those at z = 4, suggesting an even more vigorous growth in SMBH towards the high redshifts.
We caution that the C IV -based virial BH mass estimates can be under-or over-estimated due to the blueshifted component of the C IV emission line, especially for luminous quasars with potentially large C IV blueshifts.The effects can affect the conclusions obtained above, mostly in the massive and/or high-Eddington-ratio ends.For example, if we estimate the virial BH mass by adopting the calibration in Park et al. (2017), which can produce BH mass estimates similar to those with the blueshift effect corrected using the formula in Coatman et al. (2017), the resulting intrinsic BHMF shows sharp decline towards the massive end of log M BH /M ⊙ > 9.5, and the large excess at the high-Eddington-ratio end of log λ Edd > 0 is also reduced.Here, since we do not have direct measurements on the C IV blueshifts for our sample, and the calibration in Park et al. (2017) also suffers from the incompleteness of high-mass SMBHs, we adopt the results obtained with the calibration in Vestergaard & Peterson (2006) as the nominal results.Future studies on estimating the virial BH mass of our quasar samples with the Hβ and/or Mg II emission lines can be necessary to verify the determination of BHMF and ERDF at z = 4 in this work.
Target selections of z ∼ 4 quasar candidates for spectroscopic observations

Figure 1 .
Figure 1.Left) the g − r vs. r − z color-color diagram of the identified quasars at 3 < z < 4.5.Red (c1), blue (c2), magenta (c3), yellow (c4) and green (c5) shaded areas denote the c1-5 selection criteria.Black solid line shows the mean color track of the model quasars in Akiyama et al. (2018), with black squares implying its colors at z = 2, 2.5, 3, 3.5, 4, and 4.5 from left to right.Stars and dots show the identified quasars with QOP=4 (good quality) and QOP=3 (low quality), respectively.Symbols are color-coded with their spectroscopic redshift, as suggested by the color bar on the right side.Black open circles mark quasars selected by the Lyman break criteria in Ono et al. (2018) or radio detection in Yamashita et al. (2018).Right) the g − r vs. r − z color-color diagram of the identified contaminants.Orange squares, green triangles and purple inverted triangles represent the z < 3 quasars, galaxies and Galactic stars identified in the spectroscopic observations, respectively.

Figure 3 .
Figure 3. Distribution of the i-band magnitude and spectroscopic redshift, zspec, of the identified quasars.Red stars and pink dots indicate the identified 3 < z < 4.5 quasars with QOP=4 and QOP=3 quality flags, respectively.Black contours show the distributions of the luminous SDSS DR7 quasars at 3 < z < 4.5.Filled red stars and grey dots represent the 52 QOP=4 quasars and 1,462 SDSS quasars at 3.5 < z < 4.25 used in the determination of the z = 4 BHMF, respectively.Horizontal and vertical histograms plot the distributions of redshift and i-band magnitude of the quasar samples.Red open histogram, red filled histogram, pink open histogram, black open histogram and grey filled histogram represent the distributions of the 80 QOP=4 quasars, the 52 QOP=4 quasars for determining the z = 4 BHMF, the entire 151 QOP=4 plus QOP=3 quasars, the SDSS DR7 quasars at 3 < z < 4.5, and the 1,462 SDSS quasars for determining the z = 4 BHMF, respectively.Histograms of the SDSS quasars are scaled down by a factor of 45. photo

Figure 4 .
Figure 4. Comparison of the photometric redshifts of the c1 quasars with their spectroscopic redshifts.Symbols have the same meaning with Figure 3. Cyan open squares mark the quasars with |z phot − zspec| > 0.35.Black dashed line indicates the one-to-one relation.

Figure 5 .
Figure 5. Example of a fitting of the broad CIV emission line.Upper panel shows the continuum-subtracted spectrum (gray line), the best-fit Gaussian components (blue lines) and the best-fit model of the entire profile (red line) around CIV emission line.Lower panel displays residuals of the fitting by gray line.Horizontal dashed lines indicate the 1σ noise.

Figure 6 .
Figure 6.Comparison of L1350 of 80 QOP=4 quasars derived with the spectra and those derived with the broadband photometries.Symbols have the same meaning with Figure 3.The dashed line indicates the one-to-one relation, and the dotted lines delimit the area where | log L 1350−phot − log L1350−spec| < 0.3.

Figure 7 .
Figure 7. Distribution of L1350 of 80 QOP=4 quasars at 3 < z < 4.5 as a function of FWHM of CIV 1549 Å emission line.Symbols have the same meaning with Figure 3. Black contours show the same distributions for the luminous SDSS DR7 quasars at 3 < z < 4.5.Black dashed lines show a constant BH mass of 10 7 , 10 8 , 10 9 , 10 10 M⊙ from left to right.Grey shaded area suggests FWHM< 1000 km s −1 .

Figure 8 .
Figure 8. Distributions of virial BH mass (left), Eddington ratio (middle) and bolometric luminosity (right) of the quasar samples at z ∼ 4. Red and black open histograms show the distributions of the 80 QOP=4 quasars and SDSS DR7 quasars at 3 < z < 4.5, respectively, while red and grey filled histograms represent those of the 52 QOP=4 quasars and 1,462 SDSS DR7 quasars at 3.5 < z < 4.25 used in the determination of the z = 4 BHMF, respectively.Median value of each distribution is displayed with an arrow at the top.Histograms of the luminous quasars are scaled down by a factor of 65 for clarity.

Figure 9 .Figure 10 .
Figure 9. L bol vs. MBH (left) and λ Edd vs. MBH (right) distributions of the 80 QOP=4 quasars at 3 < z < 4.5.Symbols have the same meaning with Figure 3. Black contours show the same distributions of the SDSS DR7 quasars in the same redshift range.Black dashed lines in the left panel indicate a constant Eddington ratio of 1, 0.1 and 0.01 from top to bottom.Black dotted lines in the right panel denote the luminosity limits of the luminous (upper) and less-luminous (lower) quasar samples used in the determination of the z = 4 BHMF.The average 1σ uncertainties of the estimates are displayed with the gray and red errorbars for the luminous and less-luminous quasar samples, respectively.

Figure 11 .
Figure 11.Left) fraction of the less-luminous quasars among the entire c1 quasar sample.Red stars and blue dots represent the less-luminous and the entire c1 quasar samples, respectively.Grey contours show the fraction of number counts between the two samples with the color code displayed on the right.Right) effective survey area Ω(M1450, z) of the less-luminous quasar sample at z ∼ 4. Red stars and grey dots represent the less-luminous and luminous quasar samples used in the determination of the z = 4 BHMF in this work, respectively.Orange contours show the effective survey area of the less-luminous quasar sample with the color code displayed on the right.Black dashed line indicates the absolute magnitude cut of the luminous quasar sample.

Figure 12
Figure12.z = 4 broad-line AGN BHMFs (upper panels) and ERDFs (bottom panels) determined by the less-luminous (left panels), luminous (middle panels) and combined (right panels) quasar samples.Data points indicate the binned results obtained with the Vmax method, and the lines show the intrinsic results derived with the maximum likelihood method.There are four different functional forms adopted in the fitting: solid lines are for the double power-law BHMF and the log-normal ERDF, dashed lines are for the double power-law BHMF and the Schechter ERDF, dashed-dotted lines are for the Schechter BHMF and the log-normal ERDF, and dotted lines are for the Schechter BHMF and the Schechter ERDF.The best-fit results are shown by orange, grey and red lines for the less-luminous, luminous and combined quasar samples, respectively.In the right panels, pink shaded areas indicate the 1σ uncertainty of the best-fit model shown by the red solid lines.Green solid lines represent the best-fit results by adopting the log-normal ERDF with the double power-law BHMF, but excluding the mass dependence of ERDF.Red solid lines in the right panels are also plotted in the left and middle panels for comparison.

Figure 13 .
Figure 13.Comoving effective survey volume, V eff , in the λ Edd vs. MBH plane for the luminous (grey shaded area) and less-luminous quasar samples (orange contour).Symbols have the same meanings as Figure 11.

Figure 16
Figure 16.z = 4 broad-line AGN ERDFs in the two mass bins.Grey dots and orange squares show the binned ERDFs in the high mass bin of log MBH/M⊙ > 9.5 and the low mass bin of log MBH/M⊙ < 9.5, respectively.Grey and orange lines indicate the intrinsic ERDFs in the high and low mass bins, respectively.Types of lines have the same meanings as Figure 12.

Figure 17 .Figure 18
Figure 17.Observed distribution function Ψo(MBH, λ Edd ) of the best-fit double-power-law BHMF with the log-normal (left) and Schechter ERDF (right).Contours indicate the respective model distributions.Red stars and grey dots represent the less-luminous and luminous quasar samples, respectively.

Figure 19 .
Figure 19.Broad-line AGN BHMFs (left) and ERDFs (right) at z ∼ 4. Red lines show the best-fit results obtained by this work.Types of lines are the same with Figure 12.Blue and olive lines display the z = 3.75 BHMFs and ERDFs constrained by Shen & Kelly (2012) and Kelly & Shen (2013), respectively.Shaded areas indicate the 1σ uncertainty of the respective constraints.Black vertical dashed lines show the upper values where the completeness of the SDSS quasar sample drops below < 10%, given by Kelly & Shen (2013).Normalization of the ERDFs are re-scaled by integrating the corresponding BHMFs over 8 < log MBH/M⊙ < 11.

Figure 20 .
Figure 20.Cosmic evolution of the number density of broad-line AGNs at a constant BH mass (top) and Eddington ratio (bottom) at 0 < z < 6. Red dots, blue squares, purple diamonds, brown inverted triangles, light-blue open pluses, olive pluses, gray pentagons and orange stars show the results obtained in this work, Schulze & Wisotzki (2010), Schulze et al. (2015), Willott et al. (2010), Shen & Kelly (2012), Kelly & Shen (2013), Nobuta et al. (2012) and Wu et al. (2022), respectively.Green open pentagons represent the z ∼ 4 results by adopting σVM = 0.4 dex, and magenta open diamonds show those derived with the single-epoch virial BH mass calibration in Park et al. (2017).Dashed lines connect the blue squares at z = 0, purple diamonds at z = 1 − 2 and red dots at z ∼ 4 in the lowest BH mass and Eddington ratio bins; and the dotted lines connect those in the higher bins.All the ERDFs are re-scaled to match the number density of quasars in the mass range of 8 < log MBH/M⊙ < 11 before the comparison.
Figure 21.Obscuration-corrected BHMFs (top) and ERDFs (bottom) at z = 4. Black and red lines show the results with and without the correction for the fraction of obscured AGNs.Solid and dashed lines indicate the results given by the log-normal and Schechter ERDF, respectively.
shows the sharpest decline in the low-mass range of log M star /M ⊙ < 11.The resulting total BHMFs are plotted by thin lines in the left panel of Figure 22.Overall, the result obtained with the M BH − M star relation in Suh et al. (the high-z sample; 2020) is consistent with that derived with the M BH − M Bulge relation in McConnell & Ma (2013).The other two results determined with the relations in Suh et al. (the total AGN sample; 2020) and Shankar et al.

Figure 22 .
Figure 22.Left) total BHMFs (thin lines), total active BHMF (black thick line), and broad-line AGN BHMF (red thick line) at z = 4. Olive, orange, red and pink lines display the results derived with the MBH − M Bulge or MBH − Mstar relations in McConnell & Ma (2013), Suh et al. (the total AGN sample; 2020), Suh et al. (the high-z sample; 2020) and Shankar et al. (2016), respectively.The solid lines indicate the mass ranges where the MBH − M Bulge or MBH − Mstar relations are reliably constrained by the respective sample, while the dotted lines show the extrapolated ranges.Middle and right) active fraction of the broad-line AGNs (middle) and total AGNs (right) as a function of BH mass.Schulze et al. (2015) evaluate the same active fraction at z = 2 (brown dashed lines).Grey and orange shaded areas present the active fractions of the luminous and less-luminous quasars among halos at z = 4, constrained by the clustering analysis in Shen et al. (2007) and He et al. (2018), respectively.

Figure 23 .
Figure23.Luminosity functions of quasars derived by solving the continuity equation from z = 6 to 4. Results obtained with the radiative efficiency of 0.15, 0.1 and 0.05 are displayed from left to right.In each panel, red, orange and blue lines show results at z = 4, 5 and 6, respectively.Solid and dashed lines plot results obtained with the constant and time-evolving ERDF models, respectively.For the comparison, the luminosity functions of quasars at z = 4, 5 and 6 fromAkiyama et al. (2018),Niida et al. (2020) andMatsuoka et al. (2018) are plotted by red, orange and blue thin dashed-dotted lines in each panel, respectively.

Figure 24 .
Figure 24.Total BHMFs (top) and broad-line AGN BHMFs (bottom) derived by solving the continuity equation from z = 6 to 4. Results obtained with the radiative efficiency of 0.15, 0.1 and 0.05 are plotted from left to right.Lines are the same with Figure 23.For the comparison, the total BHMFs convolved from the stellar mass function at z = 4 are displayed by thin dashed-dotted lines in the top panels, and the intrinsic broad-line AGN BHMF constrained in section 4.4 is plotted by thin red dashed-dotted line in the middle panels.
at large blueshifts of > 2000 km s −1 .We derive the virial BH mass of both the luminous and less-luminous quasar samples with the calibration

Figure 25 .
Figure 25.Bivariate distributions of z ∼ 4 quasars on the L bol vs. MBH plane.Red filled stars and black dots represent the less-luminous and luminous quasar samples used to determine the z = 4 BHMF in this work.Open stars and dots show the same less-luminous and luminous quasars, respectively, but with the virial BH mass estimated with the calibration in Park et al. (2017).Black dashed lines denote constant Eddington ratios of 1, 0.1 and 0.01 from top to bottom.

Figure 26 .
Figure 26.Binned (diamonds) and intrinsic (lines) broadline AGN BHMFs (top) and ERDFs (bottom) at z = 4. Magenta color presents those with the virial BH mass recalculated with the calibration provided in Park et al. (2017).Types of lines have the same meaning with Figure 12.As a comparison, the best-fit BHMF and ERDF model using the calibration in Vestergaard & Peterson (2006) are plotted by red solid lines.
2. SMBHs, especially the massive ones, are likely to keep a high active fraction of ∼ 10% across z = 2 ∼ 4. Finally, we examine the time evolution of broad-line AGN BHMF between z = 4 and 6 through solving the He et al. continuity equation.The time evolution model of the broad-line AGN BHMF is optimized to reproduce the observed luminosity functions of quasars at z = 4, 5, and 6.We reach the following results: 1.Small radiative efficiencies of ϵ ≲ 0.1 are indicated to fully reproduce the observed luminosity functions of quasars at z = 4 and 5.

Table 2 .
Target fields of the z ∼ 4 quasars in the HSC-SSP S16A-Wide2 dataset

Table 3 .
Target fields of the z ∼ 4 quasars in the HSC-SSP S16A Udeep/Deep dataset a It adopts a different fibre configuration for the same field due to conflicts of fibre in the FoV.

Table 4 .
Identification of the z ∼ 4 quasar candidates

Table 6 .
Quasar pairs identified at z ∼ 4 in this work

Table 7 .
Continuum and emission line properties of the 80 QOP=4 quasarsNote-Table7is published in its entirety in the machine-readable format.A portion is shown here for guidance regarding its form and content.

Table 8 .
Binned broad-line AGN BHMFs at z = 4 it considers the correction for uncertainties in the virial BH mass estimates and the bolometric correction.The mass dependence of ERDF and/or the redshift evolu-tion in the BHMF and ERDF can also be flexibly taken into account.We determine the intrinsic BHMFs and

Table 10 .
Intrinsic broad-line AGN BHMFs and ERDFs determined with the combined quasar sample at z = 4