On the Connection between Supermassive Black Hole and Galaxy Growth in the Reionization Epoch

The correlation between the mass of supermassive black holes (SMBHs; $\mathcal{M}_{\rm BH}$) and their host galaxies ($\mathcal{M}_\star$) in the reionization epoch provides valuable constraints on their early growth. High-redshift quasars typically have a $\mathcal{M}_{\rm BH}$/$\mathcal{M}_\star$ ratio significantly elevated in comparison to the local value. However, the degree to which this apparent offset is driven by observational biases is unclear for the most distant quasars. To address this issue, we model the sample selection and measurement biases for a compilation of 20 quasars at $z\sim6$ with host properties based on ALMA observations. We find that the observed distribution of quasars in the $\mathcal{M}_{\rm BH} - \mathcal{M}_\star$ plane can be reproduced by assuming that the underlying SMBH population at $z\sim6$ follows the relationship in the local universe. However, a positive or even a negative evolution in $\mathcal{M}_{\rm BH}$/$\mathcal{M}_\star$ can also explain the data, depending on whether the intrinsic scatter evolves and the strength of various systematic uncertainties. To break these degeneracies, an improvement in the accuracy of mass measurements and an expansion of the current sample to lower $\mathcal{M}_{\rm BH}$ limits are needed. Furthermore, assuming a radiative efficiency of 0.1 and quasar duty cycles estimated from the active SMBH fraction, significant outliers in $\mathcal{M}_{\rm BH}$/$\mathcal{M}_\star$ tend to move toward the local relation given their instantaneous BH growth rate and star formation rate. This may provide evidence for a self-regulated SMBH-galaxy coevolution scenario that is in place at $z\sim6$, with AGN feedback being a possible driver.


INTRODUCTION
Active galactic nuclei (AGN) or quasars, powered by mass accretion onto SMBHs, produce an enormous amount of energy that has been long-speculated to have profound impacts on galaxy evolution (e.g., King & since z ∼ 2 (e.g., Jahnke et al. 2009;Schramm & Silverman 2013;Sun et al. 2015;Ding et al. 2020;Li et al. 2021). In particular, its intrinsic scatter appears similar to the local value (e.g., Ding et al. 2020;Li et al. 2021). These results suggest that a physical coupling between SMBHs and galaxies (e.g., through AGN feedback) is likely at work to keep M BH /M relatively constant. To decipher how the relationship was first established in the early universe, a key strategy would be to measure the M BH − M relation in the reionization era (z > 6) where we are able to probe the first generation of accreting SMBHs (e.g., Inayoshi et al. 2020).
Many of the z ∼ 6 quasars discovered so far are powered by extremely massive BHs with M BH ∼ 10 9 M (e.g., Fan et al. 2000;Shen et al. 2019), and are actively forming stars with star formation rates (SFR) ∼ 100 − 1000 M yr −1 (e.g., Wang et al. 2013). Their M BH /M (where M is approximated by the dynamical mass M dyn measured from gas kinematics using ALMA) appears to be significantly offset from the local value by up to 2.0 dex, suggesting that the growth of SMBHs substantially precedes their hosts (e.g., Neeleman et al. 2021). However, these quasars are biased tracers of the underlying SMBH population, since only the most luminous quasars powered by the most massive BHs can be detected in shallow surveys (e.g., Lauer et al. 2007;Volonteri & Stark 2011;Schulze & Wisotzki 2014). Lower-luminosity quasars detected in deeper surveys (e.g., the SHELLQs survey; Matsuoka et al. 2016) lie closer to the local relation thus confirming this bias (e.g., Izumi et al. 2019Izumi et al. , 2021. Moreover, the mass measurements at high redshifts suffer from significant uncertainties with possibly systematic biases. For instance, M BH of a flux-limited quasar sample might be statistically overestimated by the single-epoch virial estimator (e.g., Vestergaard & Osmer 2009; hereafter VO09) because of uncorrelated variation between AGN luminosity and broad line width (especially for Mg II and C IV). This gives rise to a luminosity-dependent virial BH mass bias (hereafter the SE bias) as detailed in Shen (2013). In addition, the hosts of z ∼ 6 quasars are found to be gas-rich (e.g., Decarli et al. 2022), thus approximating M by M dyn is likely an overestimate.
In this Letter, we model the selection and measurement biases for a z ∼ 6 quasar sample compiled in the literature in order to reveal the underlying connection between SMBH and galaxy growth in the early universe. This work assumes a cosmological model with Ω Λ = 0.7, Ω m = 0.3, and H 0 = 70 km s −1 Mpc −1 . The stellar mass and SFR are based on the Chabrier (2003) initial mass function.

SAMPLE
We adopt the z ∼ 6 quasar sample compiled in Izumi et al. (2021) while Eddington-limited accretion is assumed to estimate the mass for the remaining 24 objects. The total IR luminosity (L TIR ) is derived by fitting the 1.2 mm ALMA continuum with an optically thin graybody spectrum assuming a dust temperature of 47 K and a dust spectral emissivity index of 1.6 that have been regularly adopted for z ∼ 6 quasars, then extrapolating to the total IR (8−1000 µm) range; although the actual dust temperature could vary from source to source (e.g., Venemans et al. 2016). The SFR is derived using 3.88 × 10 −44 L TIR (Murphy et al. 2011), assuming that the cold interstellar medium is predominantly heated by star formation.
The dynamical masses M dyn of these quasars are derived through gas kinematics using the [C II] line. The standard rotating thin disk approximation is assumed for all quasars except two (given as upper limits). In this work, we only consider the 20 objects whose M BH and M dyn are derived from the Mg II line and the rotating thin disk assumption, respectively, to ensure relatively reliable mass measurements. However, we caution that the derived M dyn is highly sensitive to the assumptions made on galaxy geometry and inclination angle (e.g., Wang et al. 2013).
Ignoring the possibly small contribution of dark matter within the [C II] emitting region (e.g., Genzel et al. 2017), we estimate the M of these quasars by subtracting the expected molecular gas mass (M gas ) from their total M dyn , assuming that quasar hosts have similar gas content as star forming galaxies (e.g., Molina et al. 2021). We adopt the typical M gas /M ratio (µ gas ) of z ∼ 6 galaxies given by Tacconi et al. (2018): assumed δMS = 0. We then use this relationship to estimate the typical M gas at a given M and derive the correlation between M and M dyn , where M dyn is approximated by M + M gas . The resulting M gas and M of our quasars at their respective M dyn are shown in Figure 1 (left panel). The gas mass is distributed between ∼ 10 10 − 10 11 M , which is consistent with recent direct measurements in z ∼ 6 quasars (e.g., Decarli et al. 2022). The derived M is typically 0.2 − 0.5 dex smaller than M dyn . We also show the distribution of our quasars in the SFR -M plane in Figure 1  We adopt the local relation given by Häring & Rix (2004) (HR04) to be self-consistent with the VO09 virial estimator (see Section 6.2 in Ding et al. 2020 for the discussion on the choice of the local baseline). The local sample consists of massive ellipticals and bulgedominated S0 galaxies, thus we adopt M bulge ≈ M in the following analyses. This figure confirms that quasars at z ∼ 6 typically lie above the local relation, and the offset at a given stellar mass generally increases with quasar luminosity.

Estimating Expected Biases
As introduced in Section 1, the current z ∼ 6 quasar sample suffers from strong selection biases. Following Li et al. (2021), we perform a simple Monte Carlo sim-ulation to build a mock AGN sample that mimics the observational biases to account for such an effect in order to reveal the underlying connection between SMBHs and their host galaxies (e.g., Schulze & Wisotzki 2014;Volonteri & Reines 2016). In the following, we briefly introduce our simulation method. The details of each step and the choice of model parameters are described in Appendix A.
Our simulation starts with the galaxy stellar mass function (SMF) at z ∼ 6 given by Grazian et al. (2015) and the M BH − M relation to generate a sample of mock galaxies and SMBHs. The M BH − M relation is assumed to have a gaussian intrinsic scatter (σ µ ) with a mean that evolves as where q is the mass ratio M BH /M . Assuming that type 1 AGNs follow the same M BH −M relation as the underlying galaxy population (see Schulze &Wisotzki 2014 andLi et al. 2021), we determine the bolometric luminosity by using M BH and adopting an intrinsic Eddington ratio (λ Edd ) distribution function (ERDF) of type 1 AGNs. We use the ERDF at z = 4.75 given by Kelly & Shen (2013) who jointly constrained the intrinsic ERDF and the active BH mass function (BHMF) based on uniformly-selected SDSS quasars with the sample incompleteness being carefully corrected in a Bayesian framework. We then derive a virial M BH for each mock AGN using the VO09 virial estimator based on its true M BH , luminosity, an assumed FWHM distribution (in a lognormal form; Shen 2013), and a parameter β (0 β 1) that describes the fraction of correlated response of line width to the variation of luminosity. The value of β is poorly constrained at present. We adopt β = 0.6 in this work, while β = 1 will give rise to the SE bias (Shen 2013). The resulting virial M BH has a 0.4 dex scatter relative to the true M BH and it tends to overestimate the true M BH if L > L (and vice versa), where L is the luminosity of an AGN with a true BH mass of M BH,true , and L is the average luminosity of all AGNs at the same M BH,true (see Appendix A for details). In addition, we add a random gaussian error with a dispersion of 0.5 dex to each true M to reflect the large uncertainties of estimating M from M dyn .
In our framework, under the assumption of no evolution in the M BH − M relation (i.e., γ = 0.0, σ µ = 0.3) of the underlying SMBH population, one can estimate the expected bias (i.e., a positive ∆ log M BH relative to the local relation) caused by the sample selection function by applying the same selection criteria of observations to mock AGNs. However, it is infeasible to define a selection function as our sample is a mixture of z ∼ 6 quasars from various surveys with additional requirements of having near-IR spectroscopic and ALMA follow up to measure M BH and M dyn . Therefore, we assume a simplified scenario for which all selection biases come from the "effective" magnitude limit of different surveys, where effective means that these quasars are the relatively luminous ones selected from their parent samples for follow up observations. We simulate this selection function by producing randomly drawn samples of mock quasars that are matched to the observed M 1450 distribution (hereafter the Mock-Q sample).
The distribution of the Mock-Q sample in the M BH − M plane is plotted as blue contours in Figure 2. Their virial M BH tends to overestimate the true M BH by ∼ 0.25 dex (see Figure 5e in the Appendix). Given the magnitude limit and the large uncertainties being added to both masses, the distribution is strongly modulated compared to the originally assumed HR04 relation. In Figure 2 we also show the average virial M BH in bins of M for the Mock-Q sample as a blue solid curve. It represents the expected offset caused by selection effects and measurement uncertainties. To rephrase, any offset and large scatter in the observed M BH − M relation (red contours in Figure 2) that follows the blue curve and contours could be considered as lacking significant evolution in the mass relation of the underlying SMBH population, which appears to be true for our quasar sample. Note that the small offset between the observed quasars and the model predictions could be due to the different methods used to derive stellar masses in Grazian et al. (2015) and this work (SED fitting vs. M dyn − M gas where the latter may underestimate the total stellar mass; see Section 3.3).
We also show the impact of varying β in Figure 2. For reference, the choice of an extreme value (unlikely to be true; see Appendix A for details) for β will cause systematic shifts of the average M BH − M relation for the Mock-Q sample by about −0.2 dex and +0.4 dex for the respective case of β = 1.0 (i.e., no SE bias) and β = 0.0 (i.e., the maximum SE bias for which the BH mass tends to be overestimated by ∼ 0.65 dex). In both cases, the expected positions of z ∼ 6 quasars (i.e., the dashed and dotted curves) are offset from the actual observed ones, thus a more positively evolving (for β = 1.0) or a more negatively evolving (for β = 0.0) M BH − M relation is required to explain the observations.

Constraining the Intrinsic M BH − M Relation
While the observed offset can be reproduced by an unevolving M BH − M relation, evolutionary models cannot be ruled out. Therefore, we determine the constraints on the intrinsic evolution by generating mock AGNs with a range of γ and σ µ (assuming β = 0.6) and comparing the resulting M BH − M relation (incorporating observational biases) with the observed one to derive the likelihood (see Section 5 in Li et al. 2021 for details). The posterior distributions of γ and σ µ assuming a bounded flat prior (−2.0 < γ < 2.0, σ µ > 0.1) are shown in Figure 3. There is a strong degeneracy between γ and σ µ : either a positive evolution (i.e, M BH /M increases with redshift) with a small σ µ , or a negative evolution with a large σ µ can reproduce the large apparent offsets. The peak of the posterior distribution is slightly skewed towards a positive evolution with a small scatter, as the sample only probes overly massive quasars that are clustered at the top left of the local relation. The best-fit values based on the 16th, 50th, and the 84th percentiles of the marginalized posterior distributions are γ = 0.10 +0.80 −1.40 and σ µ = 0.38 +0.24 −0.20 , which is consistent with an unevolving M BH /M out to z ∼ 6 (e.g., Volonteri & Stark 2011;Schulze & Wisotzki 2014).
However, given various systematic uncertainties involved in mass measurements and model assumptions, it is improper to interpret the complex evolution with a single number based on a certain model. For instance, the simplified thin-disk approximation will yield an underestimated M dyn (thus M ) if the galaxy contains a dispersion-dominated component (e.g., Pensabene et al. 2020). In addition, while we intended to use the total M of these quasars in our analysis which has been shown to correlate better with M BH than using M bulge at high redshifts (e.g., Jahnke et al. 2009;Schramm & Silverman 2013;Li et al. 2021), the [C II] emission line mainly traces the inner galaxy region but not the entire galaxy (e.g., Venemans et al. 2017). These effects, together with the uncertain gas fraction and the choice of β, can induce systematic shifts of the mass measurements and affect the posterior distributions. Moreover, the bias estimates also depend on the assumed underlying distribution functions (e.g., the ERDF), which are not well-constrained at z ∼ 6 (see Section 6.3 in Li et al. 2021 for a detailed discussion). Therefore, it is clear that the current dataset is insufficient to robustly constrain the intrinsic M BH − M relation of the underlying SMBH population at z ∼ 6.

Predicted Evolution in the M BH − M Plane
Although the intrinsic evolution is unclear, overly massive quasars do exist. It is inevitable that their vigorous BH growth needs to be inhibited at some point in order to avoid unreasonably large M BH and to prevent them from further deviating from the local relation. The subsequent locations of these quasars in the M BH − M plane could be assessed by their instantaneous BH growth rate (BHGR) and SFR (e.g., Venemans et al. 2016). The BHGR can be estimated aṡ where c is the speed of light and η = 0.1 is the assumed radiative efficiency. The adopted constant η is suitable for our moderately accreting SMBHs that span 0.15 λ Edd 3.0 as expected from standard thin disk theory (Shakura & Sunyaev 1973) or slim disk theory (Abramowicz et al. 1988) for mildly super-Eddington quasars (e.g., Inayoshi et al. 2019Inayoshi et al. , 2020. We assume that these z ∼ 6 quasars can continue to form stars at their current SFR for a period (∆t). At the same time, the SMBHs keep accreting at the measured BHGR for a fraction of this time (i.e., the AGN duty cycle). The challenge is to constrain how long a quasar can sustain its high growth rate with its own feedback (e.g., Valentini et al. 2021); recent observations report a short quasar lifetime at z ∼ 6 (∼ 10 6 yr on average; Eilers et al. 2021). We estimate the M BH -dependent duty cycle from the active SMBH fraction (see the inset in Figure 4), which is derived from the ratio between the BHMF of type 1 AGNs at z ∼ 4.75 (Kelly & Shen 2013) to the total BHMF scaled from the SMF (Grazian et al. 2015) at the same redshift and assuming a M BH − M relation with γ = 0.1 and σ µ = 0.38 (Section 3.3). This method likely yields an upper limit on the BH growth time at the quasar accretion rate as most active BHs at a given M BH are of much lower luminosities. We adopt ∆t as the minimum value between the gas depletion timescale (t del = M gas /SFR; ∼ 10 − 1000 Myr) and 100 Myr (arbitrary value chosen to better visualize the evolution trend, which is shorter than most t del ). The direction of M BH and M during this period are illustrated by the dashed arrows in the top panel of Figure 4, with the caveats that our estimation is a simplification of the complex physical processes (e.g., AGN feedback, gas accretion from the environment, merger) that could happen over the next ∼ 100 Myr and the duty cycles for individual quasars are prone to significant uncertainties that are impossible to accurately constrain at present. Interestingly, the predicted evolution exhibits a flow pattern, where quasars that are significant outliers in M BH /M tend to converge to the local relation as previously seen at z < 2 (e.g., Merloni et al. 2010;Sun et al. 2015). We also show the evolution vectors (see the solid arrows) after correcting the active fraction for type 2 AGNs (see the inset in Figure 4) assuming the luminosity-dependent obscured fraction (∼ 60 − 80% at z ∼ 6) given by Vito et al. (2018) to account for the possible underestimate of the BH growth time that is not captured by the BHMF of type 1s. Still, the general trend remains. The converging pattern also holds if we adopt the peak value of γ and σ µ when estimating the total BHMF which results in lower AGN duty cycles. However, the SFR derived from the single-band ALMA photometry may be overestimated since quasars could also contribute to rest-frame ∼ 158 µm emissions (e.g., McKinney et al. 2021). Taking this effect into account will make the converging trend weaker or even disappear if the SFR is overestimated by a factor of 2. Multiband ALMA photometry are thus crucial to accurately constrain the dust temperature and subtract the quasar contamination when deriving the SFR.
In the bottom panel of Figure 4 we plot the intersection angle θ between the evolution vectors derived from the type 1+2 duty cycle and the local relation as a function of offset (color-coded by M BH ) where a decreasing trend is evident. There are 12 of 14 quasars with an offset larger than 1.0 dex have θ < 48 • , where θ = 48 • corresponds to the slope of the local relation. The median intersection angle at offset > 1.0 dex is ≈ 23 • , which is significantly smaller than that at offset < 1.0 dex (θ ≈ 54 • ). It can also be seen that at similar BH masses, θ is smaller for quasars with larger offsets, thus the decreasing trend is not driven by less massive BHs with smaller offsets and shorter duty cycles. A natural explanation of the flow pattern and the decreasing trend could be AGN feedback (e.g., Valentini et al. 2021), which sup-  Figure 4. Top: predicted evolution of the growth of SMBHs and their host galaxies over the next min(t del , 100) Myr. The dashed arrows represent the evolution vectors assuming the AGN duty cycle derived from the type 1 BHMF, while the solid arrows show the evolution vectors after correcting the duty cycle for type 2 AGNs. Two objects with upper limits in SFR are shown as black (gray) arrows. The adopted duty cycle -MBH relations are shown in the inset. The blue and orange arrows show the possible evolution pathways for overly massive quasars. The black dashed line and the shaded region represent the local HR04 relation and its scatter. Bottom: the intersection angle as a function of offset. Each quasar is color-coded by their BH masses. Two objects with upper limits in SFR are marked by blue arrows. The dashed line at θ = 48 • represents the slope of the local relation.
presses the growth of SMBHs once they deviate significantly from the local relation. The converging pattern for the most luminous and massive BHs may suggest that they have experienced rapid acccretion episodes during seeding epochs and remain being overly massive until reaching the local relation (i.e., path A in Figure 4) as shown by recent numerical simulations (e.g., Inayoshi et al. 2022). However, we cannot rule out the possibility that their progenitors are low-mass BHs moving upwards at similar stellar masses (i.e., path B in Figure 4). Such low-M BH objects are undersampled in current surveys due to the detection limit, and their vigorous BH accretion may occur rapidly in a highly obscured or/and a radiative inefficient mode which further reduces their apparent luminosity (e.g., Trebitsch et al. 2019;Davies et al. 2019). Therefore, it is essential to study the growth of low mass systems and obscured quasars that have recently been discovered at z ∼ 6 (e.g., Onoue et al. 2021).

CONCLUDING REMARKS
The z ∼ 6 quasars typically have M BH /M significantly larger than the local value. However, strong selection biases and significant measurement uncertainties severely limit the interpretation of the data. In this work, we account for these factors and demonstrate that the large apparent offsets and observed scatter could be reproduced by assuming that the underlying SMBH population at z ∼ 6 follows the local M BH − M relation ( Figure 2). However, a positive or even a negative evolution can also explain the data, depending on the evolution of the intrinsic scatter and various systematic uncertainties (Figure 3). It is thus crucial to emphasize that the evolution of the offset cannot be properly assessed without considering the scatter (see Li et al. 2021 for a similar issue at z < 1). Interestingly, quasars that are significant outliers in M BH /M tend to have evolution vectors pointing toward the local relation (Figure 4). This may provide evidence that a self-regulated SMBH-galaxy coevolution scenario is already in place at z ∼ 6, possibly driven by AGN feedback, although a robust conclusion can only be achieved with future observations that can accurately constrain the SFR (currently estimated from a single-band ALMA photometry assuming that the cold interstellar medium is mainly heated by star formation) and duty cycle for these quasars.
To break the degeneracy, expanding the current sample in both number statistics and to lower M BH limits are imperative (e.g., Habouzit et al. 2022). This is ex-pected to be achieved by the ongoing SHELLQs survey (e.g., Matsuoka et al. 2016) and the forthcoming surveys by the Vera C. Rubin Observatory and Euclid, which will offer promisingly large and less-biased quasar samples with a more uniform selection function. It is also crucial to reduce the uncertainties (especially systematic effects) in the mass measurements. The James Webb Space Telescope will enable us to measure M BH using the more reliable Hβ line, and makes it possible to probe the stellar emissions of z ∼ 6 quasars in the rest-frame optical bands thus allowing direct measurements of their M (e.g., Marshall et al. 2021). In the meantime, extending the current reverberation-mapped AGN sample to a wider parameter space is also important to validate and improve the virial M BH estimator, which will be achieved by the ongoing SDSS-V Black Hole Mapper survey. These efforts, together with a deeper understanding of the AGN accretion process (e.g., radiative efficiency, duty cycle) will allow us to better assess the connection between SMBH and galaxy growth in the reionization era. bolometric corrections to be 5.15 and 4.4, respectively (Richards et al. 2006). After these steps, we have the full knowledge of the true distribution of mock AGNs in the M − M BH − L plane. We then add realistic uncertainties to the mass terms to resemble observations. Given the large uncertainties on M dyn and M gas , we first add a random gaussian uncertainty with a standard deviation of 0.5 dex to each true M (Figure 5c). We also derive a virial M BH for each mock AGN using their true M BH , L 3000 , and an assumed Mg II FWHM distribution through Equation 1. In this step, we take into account the luminosity-dependent SE bias. The SE bias originates from uncorrelated scatter between AGN luminosity and broad line width due to both the variability of an individual quasar and the object-by-object diversity in the broad-line region (BLR) properties at a fixed true M BH (Shen 2013). We consider the following cases to represent different levels of the SE bias: • Case A: virial M BH is an unbiased estimator of the true M BH regardless of AGN luminosity. This is done by assuming that the FWHM of Mg II follows a log-normal distribution with the mean value determined by the true M BH and L 3000 for each mock AGN using Equation 1. We randomly sample the log-normal distribution with a dispersion of σ FWHM to generate FWHM for each source. The sampled FWHM is then combined with L 3000 to derive virial M BH . The dispersion σ FWHM is chosen such that the resulting scatter of virial M BH to true M BH is 0.4 dex. By doing so, the variation of luminosity (relative to the mean) at a fixed true M BH is compensated by the concordant variation in FWHM.
• Case B: Only part of the variation in luminosity can be compensated by line width. To simulate such a situation, for each mock AGN, we derive the difference (∆L) between its luminosity (L) and the mean luminosity (L) for all mock AGNs of the same true M BH . We assume that a fraction (β) of ∆L can be compensated by the line width, by using L+β∆L to determine the mean of the log-normal FWHM distribution for each source. In this case, the higher (lower)-than-the-mean luminosity can only be partly compensated by the lower (higher)-than-the-mean line width, thus the resulting virial M BH tends to overestimate (underestimate) the true M BH except at L = L.
It is currently unclear how strong β is. The non-breathing effect of Mg II (i.e., the broad line width does not respond to the continuum variability in individual quasar) suggests that β is not one (e.g., Yang et al. 2020). However, despite the lack of a BLR size (R) − L relation for individual quasars (in case of Mg II), a global R − L relation for a population of quasars spanning a broad range in BH masses and luminosities may still exist (e.g., Homayouni et al. 2020). This justifies the foundation of using the Mg II line as a single-epoch virial estimator, thus β is not likely to be zero. We adopt a relatively high response fraction (β = 0.6) as our fiducial model. This assumption yields ∆log FWHM Mg II ∝ −0.15∆L 3000 , as expected if the slope of the R − L relation for Mg II is ∼ 0.3 (e.g., Homayouni et al. 2020 With the aforementioned steps, we have generated M and virial M BH for each mock AGN with realistic uncertainties (Figure 5f). In Figure 5 we show the distributions of a mock AGN sample matching in M 1450 with the observed z ∼ 6 quasars in red. The originally assumed underlying distributions are strongly modulated by the magnitude limit and the large uncertainties on both masses. In particular, the virial M BH tends to overestimate the true M BH by ∼ 0.25 dex for this specific mock sample (Figure 5e), and the virial M BH vs. M relation is significantly offset from the HR04 relation (Figure 5f).