Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs). XVII. Black Hole Mass Distribution at z ∼ 6 Estimated via Spectral Comparison with Low-z Quasars

We report the distribution of black hole (BH) masses and Eddingont ratios estimated for a sample of 131 low luminosity quasars in the early cosmic epoch (5.6 < z < 7.0). Our work is based on the Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs) project, which has constructed a low luminosity quasar sample down to M 1450 ∼ − 21 mag, exploiting the survey data of Hyper Suprime-Cam installed on Subaru Telescope. The discovery spectra of these quasars are limited to the rest-frame wavelengths of ∼1200–1400 Å, which contain no emission lines that can be used as BH mass estimators. In order to overcome this problem, we made use of low-z counterpart spectra from the Sloan Digital Sky Survey, which are spectrally matched to the high-z spectra in overlapping wavelengths. We then combined the C iv emission line widths of the counterparts with the continuum luminosity from the SHELLQs data to estimate BH masses. The resulting BH mass distribution has a range of ∼107–10 M ⊙, with most of the quasars having BH masses ∼108 M ⊙ with sub-Eddington accretion. The present study provides not only a new insight into normal quasars in the reionization epoch, but also a new promising way to estimate BH masses of high-z quasars without near-infrared spectroscopy.


INTRODUCTION
Measuring BH masses of high-z (z ≳ 6) quasars is essential to understanding the evolution of supermassive black holes (SMBHs) and their host galaxies, as well as the mechanism of mass accretion and quasar radiation in the early universe.High-z quasar surveys have been carried out over the past two decades with deep and wide-field observations, mainly in the optical to nearinfrared (NIR) wavelengths.Now, the number of high-z quasar samples has reached more than 300 (Fan et al. 2000(Fan et al. , 2003(Fan et al. , 2004(Fan et al. , 2006;;Jiang et al. 2008Jiang et al. , 2009Jiang et al. , 2015;;Willott et al. 2005Willott et al. , 2007Willott et al. , 2009Willott et al. , 2010a,b;,b;Bañados et al. 2014Bañados et al. , 2016;;Venemans et al. 2015;Mazzucchelli et al. 2017;Reed et al. 2015Reed et al. , 2017Reed et al. , 2019;;Yang et al. 2019;Wang et al. 2017Wang et al. , 2018Wang et al. , 2019;;Mortlock et al. 2011;Venemans et al. 2013;Fan et al. 2022) including the most distant quasars known at z > 7.5 (Bañados et al. 2018;Yang et al. 2020;Wang et al. 2021).Despite the many quasars discovered in the early universe, it is unclear how such massive objects at the center of galaxies form and grow; the SMBHs in the most luminous high-z quasars have masses of M BH > 10 9 M ⊙ in a short time (< 1 Gyr) from Big Bang.For reference, the time needed for a BH with the mass of < 10 5 M ⊙ to build up a SMBH with 10 9 M ⊙ is at least 1 Gyr, assuming Eddington-limited accretion (see a recent review by Inayoshi, Visbal & Haiman 2020).No theoretical models that can fully explain this very early formation and growth have been established yet, but there are a few candidates attracting support (e.g., Volonteri 2012;Haiman 2013).The first one assumes that the progenitors are the remnants of Pop-III stars with the mass of a few 10 -100 M ⊙ , which collapse into BH seeds and then grow very efficiently to M BH ∼ 10 5 M ⊙ (e.g.Loeb & Rasio 1994;Oh & Haiman 2002;Bromm & Loeb 2003;Begelman et al. 2006;Lodato & Natarajan 2006;Inayoshi et al. 2014;Regan et al. 2014;Becerra et al. 2015;Latif et al. 2016;Volonteri & Rees 2005).The second one is the dense star cluster scenario; a dense star cluster forms when locally unstable gas flows toward the galaxy center, then the stars merge into a very massive star and finally collapse to a BH with a few thousand solar masses (e.g., Portegies Zwart et al. 2004;Omukai et al. 2008;Devecchi & Volonteri 2009).The third one is the direct collapse scenario, which starts from a supermassive star formed within the globally unstable galactic gas disks and leaves heavy BH seeds with up to ∼ one million solar masses (e.g., Begelman et al. 2006;Agarwal et al. 2012).Once formed, those seeds grow in mass via gas accretion and mergers through the hierarchical structure formation, and could produce observed high-z quasars.
In order to disentangle the above different models, we first need to obtain the distributions of BH mass and Eddington ratio to quantify SMBH assembly in the early cosmic epoch.Shen et al. (2019) (S19) have reported the properties of 50 quasars at z ≥ 5.7 based on NIR spectra collected with the Gemini telescope.The S19 sample has the BH mass range of ∼ 10 8−10 [M ⊙ ] with the average Eddington-ratio of 0.3, and no quasars have hyper massive SMBH (M BH > 10 10 M ⊙ ) or super-Eddington (λ Edd > 1) accretion.Their sample has similar UV spectral properties to a low-z luminosity-matched sample, both in the emission line profile and continuum.Yang et al. (2021)(Y21) have further studied the M BH distribution of 37 quasars at 6.30 ≤ z ≤ 7.64.The spectra were obtained with the Keck, Gemini, VLT, and Magellan telescopes.Most of their sample host massive SMBHs, which span the mass range of M BH ∼ (0.3-3.6)×10 9 M ⊙ , resulting in the predicted seeds masses larger than 10 3−4 M ⊙ on the assumption of continuous Eddington accretion since z = 30.The Fe ii/Mg ii flux ratios are comparable to low-z values and thus suggest the metal-rich environment.However, the previous high-z quasar studies are limited to the most luminous quasars with bolometric luminosity (L bol ) > 10 46 erg s −1 .We are currently carrying out a quasar survey project based on the Hyper Suprime-Cam (HSC) Subaru Strategic Program (SSP), a deep and wide-field survey with five broadband filters (g, r, i, z, y) plus narrow-band filters (Miyazaki et al. 2018;Aihara et al. 2018).The HSC-SSP survey has three layers (Wide, Deep, and UltraDeep) with different combinations of area and depth.The Wide layer covers 1200 deg 2 while the UltraDeep layer has a depth of r AB ∼ 28 mag.Based on this survey, our team, "Subaru High-z Exploration of Low-Luminosity Quasars (SHELLQs)" project has constructed a sample of 162 low-luminosity quasars at z ∼ 6-7 (Matsuoka et al. 2016(Matsuoka et al. , 2018(Matsuoka et al. , 2019a(Matsuoka et al. ,b, 2022) ) with the most distant object at z = 7.07 (Matsuoka et al. 2019a).The luminosity range of the SHELLQs sample reaches about two orders of magnitude down the previous high-z quasars sample, i.e., M 1450 < −21 mag.This low luminosity sample is expected to include low-mass SMBHs at such high redshift.On the other hand, due to their faintness, most of the SHELLQs quasars are not detected in public near-IR surveys.Onoue et al. (2019) present the measurements of Mg ii based BH mass of six relatively luminous (M 1450 ≤ −24 mag) SHELLQs quasars (three quasars also have measurements of C iv based BH masses), based on NIR follow-up spectroscopy.They found a wide black hole mass range of M BH = 10 7.6−9.3M ⊙ and the Eddington ratio (λ Edd = L bol /L Edd , where L Edd is Eddington luminosity) range of λ Edd = 0.1-1.0,but most black holes have masses of ∼ 10 9 M ⊙ and accrete with sub-Eddington ratios.It is difficult to reproduce the observed BH masses by keeping the observed sub-Eddington growth, so it implies much more active growth in the past.Thus, we may be witnessing their activity transition phase, changing from the active to quiescent mode.Izumi et al. (2018Izumi et al. ( , 2019Izumi et al. ( , 2021a,b) ,b) carried out ALMA observations of nine SHELLQs quasars and detected [C ii] (158 µm) emission line and infrared continuum.They show that the host galaxies of the SHEL-LQs sample are close to the main sequence of co-eval star-forming galaxies when they used ALMA dynamical mass as a surrogate for stellar mass, in contrast to more luminous high-z quasars.The SMBH -host dynamical mass relation in the SHELLQs sample is close to the local SMBH-stellar mass relation, which may suggest that the co-evolution is already in place at z ∼ 6.This is in clear contrast to luminous z ∼ 6 quasars, whose BH masses are above the local relation.This is the 17th paper from the SHELLQs project and is the first to present BH masses for a significant fraction of our sample.In this work, we tried to estimate their BH masses by comparing their spectra with a massive sample of low-z quasar spectra.It also provides a new way to estimate the BH mass of high-z quasars with reasonable accuracy, making it possible to discuss the evolutionary stage of the high-z low luminosity quasars.This paper is organized as follows.We describe the data and the method to get "low-z spectral counterparts" in Section 2. The main results are presented in Section 3. The growth history of our sample and some other spectral properties are discussed in Section 4. We adopt the cosmological parameters H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3, and Ω Λ = 0.7.All magnitudes in the optical and NIR bands are presented in the AB system (Oke & Gunn 1983).

DATA & ANALYSIS
We exploit a huge sample of low-z quasars and extract a "counterpart" to each of the SHELLQs quasars, with the following two assumptions.First, we assume that there are correlations between the spectral shape at λ rest ∼ 1200-1400 Å, where λ rest is rest-frame wavelength (observed by SHELLQs project), and that at the longer wavelengths which contain BH mass tracers, such as C iv, Mg ii, and Hβ emission lines.The second assumption is that there exist one or more quasars in the low-z universe that share the spectral shapes with the high-z quasars.We analyze 139 type-1 quasars at 5.66 ≤ z ≤ 7.07 with rest-ultraviolet (rest-UV) absolute magnitudes M 1450 ∼ −26 to − 21 mag (median of −24 mag), constructed by the SHELLQs project.Their spectra have been obtained by two instruments: Faint Object Camera and Spectrograph (FOCAS; Kashikawa et al. 2002) installed on the Subaru Telescope and the Optical System for Imaging and Low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS; Cepa et al. 2000) installed on the Gran Telescope CANARIAS (GTC).The two spectrographs provide spectral coverage of λ obs = 0.75 − 1.05 µm and 0.74 − 1.0 µm (where λ obs is observedframe wavelength), respectively.The spectral resolution (R) ∼ 1500 and 1200, respectively.Our discovery spectra cover only a small portion of the rest-UV wavelength at λ rest ∼ 1200 to 1400 Å, and thus do not include emission lines to derive black hole masses, e.g., C iv and Mg ii lines.
We have used SDSS Data Release fourteen (DR14) quasar catalog (Pâris et al. 2018) as a low-z quasar sample to extract the counterpart spectra.The catalog contains ∼520,000 quasars discovered over a quarter of the all-sky (∼ 10000 deg 2 ) observed through SDSS-I/II/III/IV with various selections.For example, the targets of SDSS-IV were selected by multiple selection algorithms using X-ray, optical, infrared, and radio data.X-ray sources were observed in SDSS-IV/SPIDERS (SPectroscopic IDentification of ERosita Sources).On the other hand, SDSS-IV/eBOSS selected targets based on three imaging data: SDSS, Wide-field Infrared Survey Explorer (WISE : Wright et al. 2010), and the Palomar Transient Factory (PTF; Rau et al. 2009;Law et al. 2009).Quasar candidates for SDSS-IV/eBOSS CORE sample were selected with morphological requirements and optical andWISE color cut based on XDQSOz method (Bovy et al. 2012).At higher redshift (z > 2.1), candidates were selected using their PTF photometric variability.In addition, known quasars that have low-quality spectra in SDSS-III/BOSS were re-observed in SDSS-IV.SDSS-IV/eBOSS also selected candidates from sources within 1 ′′ of a radio detection in the FIRST point source catalog (Becker et al. 1995), and also selected "time domain spectroscopic survey" targets in the g, r and i bands using the SDSS-DR9 imaging data (Ahn et al. 2012) and the multi-epoch Pan-STARRS (PS1) photometry (Kaiser et al. 2002(Kaiser et al. , 2010)).The SDSS-IV targets were observed by the eBOSS spectrographs whose resolution varies from ∼ 1300 at 0.36 µm to 2500 at 1.0 µm (Smee et al. 2013), in a series of at least three 15-min exposures.Since the SHELLQs quasars are likely a mixture of radio-loud and quiet quasars, we didn't remove radio-loud quasars from the SDSS sample when selecting counterparts.We also didn't remove BAL quasars from SDSS, since the SHELLQs sample clearly includes BAL quasars (see below).
The redshift and rest-UV magnitude ranges of the SDSS sample are 0 < z < 5 and M i [z = 2] < −20.5 mag, where M i [z = 2] represents i-band absolute magnitude at z = 2.We selected the sample with the flag "ZWARNING==0" and removed quasars that have "LOG_LYA = −999.0".100,888 quasars at 2.5 ≤ z ≤ 5.0 remain with the above criteria as a parent sample of counterparts in this work; the lower redshift cut was incorporated to include Lyα, which is the only characteristic emission line among the obtained spectra of the SHELLQs sample.Most of these low-z quasars are at z ∼ 2.5-3.0 and thus, we used their C iv lines to estimate the BH mass of SHELLQs quasars.

Extraction of "low-z counterparts"
We performed spectral fitting in λ rest ∼ 1200-1400 Å between a SHELLQs quasar and the low-z sample via χ 2 −fitting, where i represents the flux data points, f SHELLQs i and f SDSS i are the SHELLQs flux density and the SDSS flux density, respectively.A is a free parameter.The flux errors of the low-z quasars are much smaller than those of the SHELLQs quasars, and thus we considered only the errors of the latter spectra as σ i in Equation( 1).Ideally, we should select a counterpart with simillar luminosity to each of the SHELLQs quasars, but this is practically impossible since the SHELLQs quasars are significantly less luminous than the SDSS quasars at z > 2.5.We also tried to vary redshift as a free parameter, but it does not improve the results and the additional flexibility sometimes leads to catastrophic fits.Redshifts are therefore treated as the fixed values.
We extract the SDSS quasar with the most similar shape to each SHELLQs quasar.We performed χ 2 −fitting at the redward of Lyα emission line since the bluer part is severely affected by the IGM absorption for the high-z quasars.Figure 1 shows examples of the fitting results.1The spectral shapes of the low-z counterparts are indeed very similar to the features of the SHELLQs quasar spectra, not only in emission but also in absorption features; this is evident particularly in, e.g., J2216−0016, J1205−0000, J1201+0133.As we will see below, the counterpart spectra have proven to reproduce the actual spectra even at the longer wavelengths, for a limited number of objects whose NIR spectra are already available.This means that the Lyα and other spectral features in 1200 -1400 Å are indeed correlated with the spectral shape at > 1400 Å which includes the BH mass tracers.

How to measure BH masses
We estimate SMBH mass with a single epoch method based on the virial relation, using the calibration rela- , where FWHM is the full-width at half maximum of the C iv line, and λL λ is the monochromatic luminosity at 1350 Å, radiated from the accretion disc.This equation is based on the tight relation between the radius of the line-emitting region and the continuum luminosity seen in the local universe (e.g., Kaspi et al. 2005).We adopt (A, B) = (0.66, 0.53) following Vestergaard & Peterson (2006, VP06).We used the observed luminosity of each SHELLQs quasar to calculate λL λ , calculated from M 1450 assuming a power-law continuum slope of −1.5.The C iv line widths were measured from the counterpart quasars, as follows.We first subtracted the continuum emission estimated at both sides of the emission line.Next, we masked broad absorption features identified visually and removed outlying flux data points with sigma-clipping.Finally, we fitted C iv emission lines by varying the scaling factor, central wavelength, and width of single or two Gaussian profiles.The fitting wavelength range is 1450 -1650 Å.The line FWHMs were measured with those best-fit models.The uncertainties were measured with 1000 mock spectra for each counterpart generated by adding random noise to each spectral pixel based on the observed noise vectors.Figure 2 shows examples of our fitting results.We visually inspected all the fitted spectra and excluded eight quasars with poor fits, likely due to large flux errors and/or strong absorption features, from the BH mass estimates presented below.Once we estimate BH mass, the Eddington ratios are straightforward to measure.We derive the bolometric luminosity (L bol ) with L bol = 3.81 × L 1350 (Richards et al. 2006) where L 1350 is a monochromatic luminosity at 1350 Å.In reality the bolometric correction varies from one object to another, so we must keep in mind that this simple scaling from L 1350 give rise to additional uncertainty in the bolometric luminosity and Eddington ratio.We also tested obtaining 10 counterparts (with the smallest reduced-χ 2 values) and estimating BH mass as the median of the 10 measurements for each SHELLQs quasar.The comparison between this 10-counterparts case with our default case (one counterpart with the least reduced-χ 2 ) is shown in Figure 3. Overall, the two sets of measurements are consistent with each other.For simplicity and limited computing resources, we use the one-counterpart case in the following analysis and discussions.

Uncertainty of C iv based BH measurement
C iv is one of the most prominent broad BLR emission lines commonly used as BH mass estimators for z ≥ 2 quasars.However, it frequently has asymmetric profiles, broad absorptions at the blue sides of their peaks, and /or blueshift in the line centroid, which are in clear contrast to lower ionization lines such as Hβ and Mg ii (e.g., Sulentic et al. 2000;Richards et al. 2011;Brotherton et al. 2015;Gaskell 1982;Tytler & Fan 1992).Such features suggest that C iv is more severely affected than other lines by nonvirial gas motion, which would also depend on the viewing angle.The fiducial recipes of   on average consistent with the C iv estimates.However, they found a systematic offset between the C iv and Mg ii BH masses.The offset correlates with the blueshift of C iv relative to Mg ii, which suggests that the C iv line is more severely affected by a disk wind.Particularly in high luminosity quasars, C iv has been presumed to be a non-ideal line to measure BH masses, while it is the only line that can be used to estimate BH masses at high redshifts (e.g., Baskin & Laor 2005;Shen 2013;Coatman et al. 2017).On the other hand, Onoue et al. ( 2019) compared the Mg ii-and C iv-based BH masses of three high-z low-luminosity quasars and found that the two estimates show good agreement with each other.This may indicate that the C iv emitting region of lower-luminosity quasars may be less affected by gas outflows, which may be weaker.In what follows, we assume that the BH mass distributions presented below are not subject to strong systematic uncertainties while keeping all the above facts in mind.We also tested the BH mass calibration that takes into account and corrects for the effect of C iv blueshift component, presented by Coatman et al. (2017).This calibration was empirically established from 230 high-luminosity (L bol ∼ 10 45.5−48.5 erg s −1 ) quasars at 1.5 < z < 4.0, which covers both the hydrogen Balmer emission lines and the C iv emission lines.We found little change in the estimated BH masses, and thus confirmed that the conclusions of this paper remain unchanged with this alternative calibration.

Consistency of our BH mass estimates with direct measurements
Here, we compare our BH mass estimates with those based on the direct spectral measurements provided by Onoue et al. (2019, O19).O19 collected NIR spectra of six SHELLQs quasars at 6.1 < z < 6.7 and with M 1450 ≤ −24 mag, which is a practical limit to obtain reasonable emission line measurements with, say, less than 10 hours with 8-m class telescope on the ground.Specifically, O19 observed the six quasars with Very Large Telescope/X-Shooter and Gemini-N/GNIRS, with the spectral resolution of R ∼ 5000 -7000 and R ∼ 500 -800, respectively.The resultant Mg ii (λ2798) based BH masses are in the range of M BH = 10 7.6−9.3M ⊙ , with the Eddington ratio of L bol /L Edd = 0.16 − 1.1, while most of the quasars have M BH ∼ 10 9 M ⊙ with sub-Eddington accretion.Three quasars in the O19 sample also have M BH estimates based on C iv emission lines.In addition, nine more quasars from the SHELLQs sample has the observed NIR spectra to measure BH masses (M.Onoue et al., in prep.). Figure 4 shows comparisons of the continuum luminosities, line widths (FWHM), and BH masses between the direct spectral measurements (horizontal axis) and the present work based on the low-z counterparts (vertical axis).In addition to the C iv-based measurements, we include Mg ii-based measurements in this comparison; the counterparts for the latter have been newly selected from SDSS in a limited redshift range (2.1 < z < 2.5) where both Lyα and Mg ii are covered by the SDSS spectroscopy, and the BH masses have been estimated from their Mg ii in exactly the same way as we did based on C iv.These comparisons suggest that our method of estimating M BH via low-z counterparts works well with reasonable accuracy.While there are a modest amount of scatters in the plotted quantities, the BH masses from the two methods are consistent with each other within ∼ 0.3 dex, which is comparable to the typical systematic uncertainty expected in the single-epoch mass estimates.Figure 5 represent the spectral comparisons of the three quasars and their counterparts with C iv-based M BH (i.e., those corresponding to the dots in Figure 4).While there is a good overall agreement between the two sets of spectra, we observed offsets in C iv peak wavelengths of J1152+0055 and J2239+0207 from the expected positions.This may be caused by uncertainty in redshift of the SHELLQs quasars, which have been determined from the Lyα line affected severely by the IGM absorption.
Since [C ii] redshifts are available for the above three quasars from SHELLQs ALMA observations (Izumi et al. 2018(Izumi et al. , 2019)), we tested to use those more accurate redshift to convert observed spectra into the rest-frame and investigate how the redshift uncertainly affects the accuracy of this counterpart method.The counterparts and thus BH mass estimates of two out of the three quasars have changed from the case with Lyα redshift, as shown in Figure 6 (see Figure 15 in Appendix for the updated spectral comparison similar to Figure 5).We observe almost perfect agreement between the two sets of measurements in this case, demonstrating the potential power of the present method.However, counterparts with Lya redshifts (the only redshift estimates available for the majority of the SHELLQs quasars) still provide reasonable BH mass estimates with ∼ 0.2 dex difference from the [C ii] redshift cases, and such small difference does not affect the conclusion of this paper.Having said that, we stress that future follow-up observations of those quasars to measure more accurate redshifts, with [C ii], Mg ii and/or other emission lines, will significantly enhance the usefulness of this unique high-z quasar sample for various topics.
As an additional test of the present method, we randomly selected 100 quasars from z ∼ 2-3 SDSS sample and cutout their spectra in the range of λ rest =1200 -1400 Å, and then obtained their counterparts from the remaining SDSS quasars in exactly the same way as we did for the SHELLQs quasars.Figure 7 compares BH masses based on direct measurements from C iv lines contained in the original 100 spectra and those based on the counterpart method.We found a broad agreement between the two measurements, with the systematic offset and standard deviation of 0.09 dex and 0.23 dex, respectively; these results show a good correlation even on the high-mass side.
We further tested our method with luminous quasars taken from the XQ-100 sample.This sample was constructed based on the European Southern Observatory Large Programme (program code 189.A-0424) "Quasars and their absorption lines: a legacy survey of the highredshift Universe with VLT/X-shooter".The XQ-100 survey produced high-resolution (R ∼4000-7000) spectra of 100 quasars at redshift ≃ 3.5-4.5 covering Lyα and C iv emission lines (López et al. 2016).The median absolute magnitude is M 1450 = -29.6 ± 0.017 and the median bolometric luminosity is L bol = (8.0 ± 0.13) × 10 47 ergs −1 .We measured their BH mass using C iv lines with the VP06 calibration, and calculated their uncertainties based on the MCMC method with 100 mock spectra.At the same time, we estimated thier BH mass with the counterpart method using the SDSS spectra.The comparison between the two measurements is shown   in Figure 8.We found a reasonable agreement again, with the systematic offset and standard deviation of 0.13 dex and 0.22 dex, respectively.Overall, we found that our new method provides a reasonable accuracy of BH mass estimates for high-z quasars, whose near-IR spec-   The estimated BH masses and the bolometric luminosities span a wide range of M BH ∼ 10 7−10 M ⊙ and ∼ 10 45.0−47.0erg/s, with the median values of ∼ 10 8.6 M ⊙ and ∼ 10 46.1 erg/s, respectively.For comparison, we also plot the quasar sample at lower redshift, SDSS DR14 quasars with masses estimated using either C iv, Mg ii or Hβ depending on redshift in Rakshit et al. (2020).The BH masses of the luminous S19 sample at z ∼ 6 are based on either C iv or Mg ii and populate the range between 10 8−10 M ⊙ except for one quasar at > 10 10.5 M ⊙ with poor spectral fit.Our BH mass and Eddington ratio are lower on average by ∼ 1 dex and ∼ 0.15, respectively than those in the S19 sample.On the other hand, Willott et al. (2010a) (W10) reported Mg ii-based BH masses of nine high-z quasars with relatively low luminosity (M 1450 < −24.3 mag).They have M BH ∼ 10 8−10 M ⊙ with sub-Eddington to Eddington accretion.The median BH masses of W10 are not so different from our BH masses, but the median Eddington ratio is about seven times higher than our sample.Sub-Eddington accreted quasars dominate our sample, representing the unprecedented depth of our survey.Figure 10 presents a histogram of Eddington ratios in the three sample, clearly showing the dominance of sub-Eddington accretors in the SHELLQs quasars.This is consist with the prediction from a semi-analytic models of galaxy formation, reported by Shirakata et al. (2019).They calculated the Eddington ratio distributions of SMBHs with varying sample selection at 0 < z < 6, and found that the shallower luminosity cuts tend to miss AGNs with lower Eddington ratios, and that such a trend is more critical at higher redshit.
We found six quasars marginally exceeding the Eddington limit.J0859+0022 has been observed by O19, who reported the Eddington ratio of 1.1 +0.5 −0.3 .Figure 16 in Appendix shows the spectrum of one of the SHELLQs quasars with the highest estimated Eddington ratios.Their spectra tend to have high equivalents widths of Lyα and C iv, the latter being observed in the counterparts.Similarly, lower-mass (M BH < 10 8 M ⊙ ) quasars in our sample tend to have high EW of the emission lines.Figure 17 shows the spectrum of our lowest mass quasar along with that of its counterpart.Its estimated BH mass is 10 (7.0±0.2) M ⊙ and the Eddington ratio is 0.8.The high EW of the emission lines in these quasars may indicate gas-rich environments around the nuclei.Given the high Eddington ratios and/or low BH masses, they may represent the early stage of the assembly, accompanied by abundant inflowing gas. Figure 9 suggests that there are few hyper massive SMBHs ( > 10 10 M ⊙ ) in the present and other samples.This is partly due to the observing difficulties since broad emission lines become extremely broad and weak compared to the continuum in such objects.Alternatively, we could think of a possibility that the AGN timescale at the high-mass end is so short that we cannot find them in our surveys over limited portions of the sky.It is also possible that the BH mass limit is actually determined by the host galaxy masses.

SMBH growth history
The growth of a SMBH is exponential if it keeps a constant Eddington ratio.The timescale for a seed black hole with mass M seed to reach a given M BH is where τ is the e-folding timescale,

Gyr.
The radiation efficiency η is the factor describing how efficiently the accreting mass is converted to radiation.Bañados et al. (2018) reported the possible growth history of three luminous quasars at z = 6.3 − 7.5.They found 10 3 < M seed < 10 4 M ⊙ at z = 30, assuming that the quasars had been accreting at the Eddington limit with a radiative efficiency of 10%.This suggests that SMBHs in the early universe have large initial masses and/or sustained Eddington limit (or episodically super-Eddington) accretion.The most distant quasar currently known, J0313-1806 at z = 7.642 (Wang et al. 2021), has the most massive black hole (M BH ∼ 10 10 M ⊙ ) at z > 7, and thus poses a powerful constraint on the seed black hole mass.The estimated seed mass is M seed ∼ 10 4−5 M ⊙ with the assumptions of Eddington-limited accretion with a 10% radiative efficiency and a duty cycle of unity.The above results are contrary to the scenario assuming the seeds of Population III star remnants.
Here, we trace back in time the BH masses of the 131 SHELLQs quasars we determined in this work.The formation redshift is assumed to be z = 30 when the first stars and galaxies were thought to have formed (Bromm & Larson 2004;Bromm & Yoshida 2011).The radiation efficiency of η = 0.1 (i.e. a standard thin accretion disk; Shakura & Sunyaev 1976) is assumed.We consider two cases.The first we assume is that the SHELLQs quasars had grown at the Eddington limit,  and the second is that they had grown with a constant Eddington ratio estimated at their observed redshift.Figure 11 displays the first case; most of the seeds have the mass range of putative Pop-III star remnants (M seed = 10 − 100 M ⊙ ) already at z ∼ 20, in contrast to the cases of previously-studied luminous quasars (e.g.Mortlock et al. 2011;Wu et al. 2015;Yang et al. 2020;Wang et al. 2020).This is primarily because most of our sample have M BH < 10 9 M ⊙ , and it is possible for Pop-III seeds to create such SMBH masses without growth beyond the Eddington limit.At the same time, some of our sample have masses exceeding the range of the remnant Pop-III stars, in particular those at the highest redshift.Such objects may be formed from other seed populations such as dense star clusters or direct collapse black holes, or have grown with super-Eddington accretion.
Figure 12 displays the second case, the growth path of our sample with perhaps the more realistic assumption on the Eddington ratios.Most seeds exceed the range allowed for Pop-III remnants, which indicates either that we have to assume other seed populations or that they have grown with higher Eddington ratios than the values observed at z ∼ 6. Figure 12 may suggest that the distribution of BH mass and the growth history of the high-z low luminosity quasars are divided into the following two phases.The majority of the sample are indicated to have massive seeds even above those predicted by the direct collapse scenario with the assumption of constant observed Eddington ratios.Hence, they are likely to have grown with higher Eddington ratios in the past, and have switched to the less active phase by z ∼ 6.The others could be in the quiet stage with sub-Eddingon accretion from the seeds to the observed epochs.They may switch to the active phase if they eventually experience major evolutionary events such as galaxy mergers.Sixteen SHELLQs quasars have been observed at the sub-mm wavelength with ALMA, one of which has a companion separated by 15 kpc (Izumi et al. in preparation).This quasar may possibly switch to active mode in a relatively short timescale through a galaxy merger.

Spectral properties compared with luminosity-matched quasars
In this section, we compare the spectral properties of our sample with the local luminosity-matched sample, in order to test whether the relation between luminosity and spectral shape has redshift dependence.
Many previous studies have addressed the redshift dependence of quasar spectral properties.It is generally thought that quasars have little redshift evolution in their spectral properties (e.g., Mortlock et al. 2011;Shen et al. 2019).Fine et al. (2008) and Shen et al. (2008) found that the line widths of the virial BH mass estimators do not depend strongly on luminosity or redshift, especially for luminous quasars in their sample.On the other hand, Willott et al. (2010a) suggested that such results are caused by a wide range of Eddington ratios in the sample used for the analysis.Instead, they found a strong correlation between luminosity and Mg ii FWHM (and thus M BH ) of the high-z quasars from SDSS and CFHQS, which is presumably due to the fact that most of them are accreting at close to the Eddington limit.W10 also report that the Eddington ratio distribution has a peak at 1.07, which is considerably higher than that of the luminosity-matched sample at z = 2 (λ Edd peaks at 0.37), suggesting that typical quasars at the higher redshift have higher level of accretion activity.
We now proceed to check the rest-frame UV spectral properties of the SHELLQs sample, which represents the largest sample of z ≥ 6 quasars at the faint end.As the first step, we created a low-z control sample matched in continuum luminosity at the rest-frame 1350 Å to our high-z sample.As is well known, quasar luminosity governs many emission line properties, thus it is crucial to match in luminosities.For instance, in high-ionization lines such as C iv, EW decreases with luminosity (i.e., the Baldwin effect, Baldwin 1997), and other line profiles (e.g., velocity shift and asymmetry) also change with luminosity in a systematic manner (e.g., Richards et al. 2002).We randomly selected 50 control quasars from the SDSS DR14 quasar catalog within the luminosity difference of ∆logL = 0.2 for each of the SHELLQs quasars.The control quasars have redshifts z ∼ 1.5−5.0;thus, their spectra cover most of the major rest-frame UV lines from Lyα to Mg ii.As the second step, we created the mean composite spectra from each of the SHELLQs sample, counterpart sample, and control sample.
Figure 13 displays the composite spectra of the three samples.The Lyα of the SHELLQs and counterpart  composite appears to be strong compared to the control, while we see no significant difference in the C iv line between the counterpart and control.We measured EW and FWHM of these lines as follows.We subtracted the continuum fluxes estimated at the both sides of the lines, then fitted the residuals with two or three Gaussian profiles.The SHELLQs composite lacks the blue side of their Lyα line due to the IGM absorption, so we obtained the best-fit model with only the red side of the peak.The associated errors were estimated with 100 mock spectra created via the Markov Chain Monte Carlo method using the observed and propagated noise array.The results of these measurements are shown in Table 1.FWHM of Lyα is apparently smaller than others.On the other hand, Lyα EWs of the SHELLQs and counterpart are twice as high as the control.Since the SHELLQs and the control samples are matched in luminosity, this difference does not reflect the well-known Baldwin effect.However, interestingly, EWs of C iv show the opposite tend; the counterpart has a slightly smaller EW than the control.We suggest that the different behavior of Lyα and C iv may reflect the difference in physical condition of the emitting gas, such as the density and metalicity.The counterpart has smaller line widths than the control in both Lyα and C iv lines, but the difference is within 1sigma uncertainty.
Figure 14 compares the M BH and λ Edd distributions of the SHELLQs and the control sample.The p-values of the K-S test in the BH mass and Eddington ratio distribution are ∼ 0.008 and ∼ 0.0007, respectively, which may suggest that the basic properties of mass accretion are different between the two samples.The former sample has higher BH masses by 0.1 dex with lower Eddington ratios by 0.1 than the latter.At the face value, it may provide a hint of higher SMBH activity in the SHELLQs quasars at z ∼ 6 than the control sample, despite the fact that the quasar activity generally peaks around 1 < z < 3. Some previous studies found more prominent differences between high-z and low-z quasars (e.g.Willott et al. 2010a;Shen et al. 2019); the discrepancy from the present work may be attributed to the fact that we are studying unprecedentedly lowluminosity sample at such high redshifts.Having said that, we stress that the present comparisons are subject to strong systematic biases, such as those caused by sample definition, completeness, measurements errors, and assumptions used in the calculations.For robust detection of a 0.1-dex level difference in Eddington ratio distributions for example, one needs to resolve the complexity of those biases and correct for their combined effect, which would be a subject of future projects.

Final remarks
This paper presents a novel method to measure the distributions of BH mass and Eddington ratio among high-z quasars, applicable even to quasars with very low luminosities, such as those established in the SHEL-LQs program.Our obvious next step is to confirm the estimated properties with direct measurements of a part of the present objects, through future observations with NIR spectrographs either from the ground or using James Webb Space Telescope.Once confirmed with follow-up NIR spectroscopy, the present sample and measurements provide new insight into the early quasar evolution in the reionization epoch.
While we identified candidates of very low-mass SMBHs in our sample, the HSC imaging may have detected more similar objects.The SHELLQs spectroscopy program prioritized candidates with detection in both the HSC z -and y-bands, which contain the redshifted Lyα and continuum, respectively.However, the identified low-mass quasars (e.g., J0859+0022) tend to have a very strong and narrow Lyα and faint continuum, and thus are often undetected in the y-band.We consider them as the high-z analogs of narrow-line Seyfert 1 found in the local universe.There may be a significant number of similar objects in the remaining HSC candidates awaiting spectroscopy, which may populate the low-mass end of the M BH distribution.
This research is based on data collected at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.We are honored and grateful for the opportunity of observing the universe from Maunakea, which has cultural, historical, and natural significance in Hawaii.The data analysis was in part carried out on the open-use data analysis computer system at the Astronomy Data Center of NAOJ.This work is also based on observations made with the

Figure 1 .
Figure 1.Example of the spectral fitting between SHELLQs quasars (gray solid line) and their low-z counterparts from SDSS (red solid line).The dashed lines mark the expected positions of Lyα, N v, and C iv emission lines.
BH mass estimates summarized inShen et al. (2011) are known to have substantial systematic uncertainty of ∼ 0.4 dex, estimated from the differences among different broad-line estimators, such as C iv and Mg ii.This systematic uncertainty is often larger than the measurement uncertainty of spectral fits.In contrast, Mg ii line widths are well correlated with those of H β, and these two low-ionization line estimators usually give consistent virial masses ( e.g.,McLure & Dunlop 2004;Salviander et al. 2007).Shen et al. (2008)  compiled BH masses of ∼ 60, 000 quasars in the redshift range 0.1 ≤ z ≤ 4.5 based on the Hβ, Mg ii, and C iv emission lines.Within their sample, the Mg ii and Hβ based BH masses are

Figure 2 .
Figure 2. The best-fit models of continuum-subtracted spectra around the C iv lines.Black solid lines show the observed flux data, and gray shaded areas represent the masked region.The best-fit Gaussian profiles are in red, and FWHMs are in green.Black dashed lines represent the theoretical peak wavelength of 1549 Å.

Figure 3 .
Figure 3.Comparison of BH masses measured from one counterpart that has the smallest reduced-χ 2 (in the vertical axis) with those measured from ten counterparts (in the horizontal axis).The horizontal error bars represent 1 sigma scatter of the ten counterparts.

Figure 4 .
Figure 4. Comparisons between our estimates (vertical axis) and actual measurements (horizontal axis; Onoue et al. 2019, M. Onoue et al., in prep.) of the continuum luminosity (top), line FWHM (middle), and BH mass (bottom).The dots and triangles represnt the C iv-based and Mg iibased estimates, respectively.

Figure 5 .
Figure 5. Spectral comparison of three quasars with C ivbased BH mass estimates (corresponding to the dots in Figure 4).Gray solid lines show the NIR spectra obtained by O19.Black solid lines show the low-z counterpart spectra.Red and brown lines mark the peak wavelengths estimated by our Gaussian fitting and theoretical wavelength of the C iv lines, respectively.Yellow windows provide the extended plots around C iv (at 1500 -1600 Å).The two sets of spectra (SHELLQs vs. counterparts) are in reasonable agreement with each other, even outside of the spectral coverage (> 1400 Å) used for the fitting.

Figure 6 .
Figure 6.Same as Figure 4 (bottom), but for C iv-based BH mass estimates with Lyα redshifts (open circles) and with [C ii] redshifts (dots).

Figure 7 .
Figure 7.Comparison of BH masses between direct measurements and counterpart estimates for a randomly selected sample of 100 SDSS quasars at z ∼ 2-3.
3. BH MASS DISTRIBUTION OF HIGH-Z LOW-LUMINOSITY QUASARS We show the relationship between the BH masses and bolometric luminosity of the SHELLQs sample in Figure 9. Single epoch estimates of BH masses for highz quasars, based on calibrations established at lower redshift, are known to have large systematic uncer-tainties, whose typical value of 0.5 dex is indicated at the lower right of the figure.Our high-z and other samples positively correlate in this parameter space.

Figure 9 .
Figure 9.The SMBH mass-luminosity plane of high-z (z > 5.6) quasars.Blue dots represent our sample, whose BH masses were measured with the counterpart method.Blue open circles represent those with possible BAL features in the counterparts, with the SDSS measured BALnicity index BI > 0.0 (Pâris et al. 2018).The green and brown dots represent the samples from Shen et al. (2019) and Willott et al. (2010a), respectively.The contours show the distribution of SDSS DR14 quasars at lower redshifts.The diagonal lines show the constant Eddington ratios of L bol /L Edd = 1, 0.1, 0.01 from top left to bottom right.The typical systematic uncertainty of the MBH measurements (0.5 dex) is shown with the error bar at the lower right.

Figure 11 .
Figure 11.An estimated growth history of the SHELLQs quasars.The horizontal axis gives redshift (bottom) and time since the Big Bang (top).Solid lines show the cases where SMBHs have grown at the Eddington limit λ Edd = 1, with the dots representing the estimated BH masses at the observed redshift.The shaded regions correspond to the predicted typical mass ranges for the three scenarios, i.e., Pop-III remnants (MBH ≤ 10 3 M⊙; green), dense star clusters (MBH ∼ 10 3−4 M⊙; purple) and direct collapse BHs (MBH ∼ 10 4−6 M⊙; orange).

Figure 12 .
Figure 12.Same as Figure 11, but for the case of the Eddington ratios being fixed to the observed values.

Figure 14 .
Figure 14.Distributions of BH mass (left) and Eddington ratios (right) of the SHELLQs (blue) and the control (gray) quasars.
Figure 15.Same as Figure 5 but for the case in which the SDSS counterparts are obtained with [C ii] redshift.

Figure 17 .
Figure 17.Same as Figure 1, but for the lowest mass SHELLQs quasars with logMBH ∼ 10 7.0±0.2M⊙ and with a high Eddington ratio.

Table 1 .
Spectral measurements from the composites