THE IMPACT OF THE UNCERTAINTY IN SINGLE-EPOCH VIRIAL BLACK HOLE MASS ESTIMATES ON THE OBSERVED EVOLUTION OF THE BLACK HOLE–BULGE SCALING RELATIONS

The redshift evolution of the local scaling relations between galaxy bulge properties and the mass of the central supermassive black hole (SMBH) has important clues to the establishment of these tight relations across cosmic time. There is currently a huge effort in measuring this evolution, either in terms of the BH mass–bulge velocity dispersion (M_•–σ) relation, or the BH mass–bulge stellar mass/luminosity (M_•–M_bulge, M_•–L_bulge) relation (e.g., Shields et al. 2003; Peng et al. 2006a, 2006b; Salviander et al. 2007; Shen et al. 2008b; Woo et al. 2006, 2008; Treu et al. 2007; Jahnke et al. 2009; Greene et al. 2010; Merloni et al. 2010; Decarli et al. 2010; Bennert et al. 2010). With a few exceptions, most of these studies report a strong positive evolution in BH mass for fixed bulge properties, which reaches as high as ∼0.6 dex offset in BH mass from the local relation at the highest redshift probed (z ∼ 2).

This observed strong evolution is difficult to understand in two aspects. For the M_•–σ relation, numerical simulations based on self-regulated BH growth (e.g., Robertson et al. 2006) and other theoretical arguments (e.g., Shankar et al. 2009) favor a mild evolution in the normalization of the M_•–σ relation, in conflict with the otherwise claimed strong evolution. For the M_•–M_bulge or M_•–L_bulge relation, many of the observed hosts at earlier epoches are already bulge-dominated;² hence, unless their bulges continue to grow a considerable amount over time, it is difficult to understand how these systems would have migrated to the local relation.

An alternative explanation for this observed strong evolution is that it is caused, at least partly, by the systematics involved in deriving both host properties and BH masses. All these studies mentioned above focus on broad line quasar samples, with host properties measured from the quasar-light-subtracted images or spectra. The systematics in subtracting the quasar light as well as in converting observables (such as photometric colors or spectral properties) to physical quantities (such as bulge stellar mass) is a potential contamination to the final results. More importantly, the samples are predominately selected in quasar luminosity and hence are not unbiased samples for studies of the BH scaling relations.

One statistical bias resulting from using quasar-luminosity-selected samples was discussed in detail in Lauer et al. (2007). Because there is cosmic scatter in the BH scaling relations and because the galaxy luminosity (or stellar mass, velocity dispersion, etc.) function is bottom heavy, there are more smaller galaxies hosting the same mass BHs than the more massive galaxies. Hence, when selecting samples based on quasar luminosities, low-mass BHs are underrepresented, leading to an apparent excess in BH mass at fixed host properties. If the cosmic scatter in BH scaling relations increases with redshift and reaches ≳0.5 dex, it can in principle account for all the excess in BH mass observed for the most luminous quasars (e.g., Lauer et al. 2007; Merloni et al. 2010).

A second statistical bias, resulting from the uncertainty in the BH mass estimates, was pointed out in Shen et al. (2008b; also see Kelly et al. 2009a). In all these studies, the BH masses are estimated using the so-called virial method based on single-epoch spectra. In this method, one assumes that the broad line region (BLR) is virialized, and the BH mass can be estimated as M_vir ≈ G⁻¹RV², where R is the BLR radius and V is the virial velocity; one further estimates the BLR radius using a correlation between R and the continuum luminosity L, i.e., $R\propto L^{C_1}$, found in local reverberation mapping (RM) samples (e.g., Kaspi et al. 2005; Bentz et al. 2006, 2009), and estimates the virial velocity using the width of the broad lines. In this way, one can estimate a virial BH mass using single-epoch spectra: log M_vir = log R + 2log(FWHM) + const = C₁log L + 2log(FWHM) + const. The coefficients are calibrated empirically using RM samples or inter-calibrations between various lines (e.g., McLure & Jarvis 2002; Vestergaard & Peterson 2006; McGill et al. 2008).

Even though the virial method is widely used, these BH mass estimates based on several lines (usually Hα, Hβ, Mg ii, and C iv) have a non-negligible uncertainty σ_vir ≳ 0.4 dex, when compared to RM masses or masses derived from the M_•–σ relation (e.g., Vestergaard & Peterson 2006; Onken et al. 2004). This uncertainty must come from the imperfectness of using luminosity and line width as proxies for the BLR radius and virial velocity, i.e., there are substantial uncorrelated variations σ_L in luminosity and variations σ_FWHM in line width, which together contribute to the overall uncertainty σ_vir, i.e.,

$$\begin{equation} \sigma _{\rm vir}=\sqrt{(C_1\sigma _{L})^2 + (2\sigma _{\rm FWHM})^2}{.} \end{equation} \tag{ 1 }$$

One can then imagine that a statistical bias will arise if the underlying active black hole mass function (BHMF) is bottom heavy. In particular, the variations (uncompensated by the variations in FWHM) of luminosity at fixed true BH mass, σ_L, will scatter more lighter BHs into a luminosity bin than heavier BHs, and bias the mean BH mass in that bin. This is the Malmquist-type bias (or Eddington bias) emphasized in Shen et al. (2008a, Section 4.4), which has received little attention in the studies on the evolution of BH scaling relations to date.

In this paper, we examine the impact of this mass bias on the observed evolution in BH scaling relations with realistic models for the underlying true BHMF and quasar luminosity function (LF). In Section 2, we review the statistical mass bias discussed in Shen et al. (2008a) and demonstrate its effects with simple models; in Section 3, we consider more realistic intrinsic BHMF and quasar LF, and estimate the mass bias as functions of luminosity and redshift; we discuss the impact of this mass bias and conclude in Section 4. Throughout the paper, we use cosmological parameters Ω₀ = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s^-1 Mpc⁻¹. Luminosities are in units of erg s^-1 and BH masses are in units of M_☉. Unless otherwise specified, "luminosity" refers to the bolometric luminosity, and we are only concerned with the active BH population.

Denoting λ ≡ log L, where L is quasar luminosity (bolometric or in a specific band), m ≡ log M_• the true BH mass, and m_e ≡ log M_vir the virial BH mass, the probability distribution of virial mass m_e given true mass m is

$$\begin{equation} p_0(m_e|m)=\big(2\pi \sigma _{\rm vir}^2\big)^{-1/2}\exp \left[-\frac{(m_e-m)^2}{2\sigma _{\rm vir}^2}\right]{,} \end{equation} \tag{ 2 }$$

where σ_vir is the uncertainty of the virial BH mass.

The probability distribution of true BH mass m given virial mass m_e is

$$\begin{equation} p_1(m|m_e)=p_0(m_e|m)\Psi _M(m)\left(\int p_0(m_e|m)\Psi _M(m)dm\right)^{-1}{,} \end{equation} \tag{ 3 }$$

where Ψ_M(m) ≡ dn/dm is the true BHMF. If the sample is selected irrespective of luminosity (i.e., no flux limit) and so we can see all BHs with arbitrary m, then for a power-law true BHMF $\Psi _M(m)\propto 10^{m\gamma _M}$ and at fixed m_e, the statistical bias between the expectation value 〈m〉 and virial mass m_e is simply

$$\begin{equation} \Delta \log M_{\bullet }= m_e - \langle m \rangle ={-}\ln (10)\gamma _{M}\sigma _{\rm vir}^2\,{.} \end{equation} \tag{ 4 }$$

In reality the sample is usually selected in quasar luminosity. If the distribution of luminosity at a given true BH mass for broad line quasars is g(λ|m) where ∫g(λ|m)dλ = 1, the expectation value of m at fixed luminosity λ is

$$\begin{equation} \langle m\rangle _{\lambda }=\frac{\int\! mg(\lambda |m)\Psi _M(m)\,dm}{\int\! g(\lambda |m)\Psi _M(m)\,dm}{,} \end{equation} \tag{ 5 }$$

and the luminosity-weighted average of m in the quasar sample is

$$\begin{equation} \langle m\rangle = \frac{\int _{\lambda _1}^{\lambda _2}\Psi _L(\lambda) \langle m\rangle _{\lambda }d\lambda }{\int _{\lambda _1}^{\lambda _2}\Psi _L(\lambda)d\lambda }{,} \end{equation} \tag{ 6 }$$

where λ₁ and λ₂ are the luminosity limits in the sample, and the quasar LF Ψ_L ≡ dn/dλ is

$$\begin{equation} \Psi _{L}(\lambda)=\int \ \Psi _{M}(m)g(\lambda |m)\,dm{.} \end{equation} \tag{ 7 }$$

On the other hand, the virial BH masses depend on luminosities, and the sample averaged virial mass is

$$\begin{equation} \langle m_e\rangle =\frac{\int _{\lambda _1}^{\lambda _2}\!\Psi _L(\lambda) \frac{\int m_ef(m_e|\lambda)\,dm_e}{\int f(m_e|\lambda)\,dm_e}\,d\lambda }{\int _{\lambda _1}^{\lambda _2}\Psi _L(\lambda)\,d\lambda }{,} \end{equation} \tag{ 8 }$$

where f(m_e|λ) is the probability distribution of m_e given luminosity λ and ∫f(m_e|λ) dm_e = 1. In the absence of a theoretical model for f(m_e|λ), we adopt here an empirical recipe for f(m_e|λ): recall that the virial BH mass is expressed in terms of luminosity and FWHM. Neglecting the scatter from converting bolometric luminosity to continuum luminosity, which is typically ∼0.1 dex (e.g., Richards et al. 2006a), we have

$$\begin{equation} m_e=C_1\lambda + 2\log ({\rm FWHM}) + C_2{.} \end{equation} \tag{ 9 }$$

Here, C₁ ≈ 0.5–0.7 is the slope in the measured luminosity–radius relation for BLRs and C₂ is calibrated empirically (e.g., Kaspi et al. 2005, 2007; Bentz et al. 2006, 2009; McLure & Jarvis 2002; Vestergaard & Peterson 2006). For statistical quasar samples, the distribution of FWHMs does not depend much on luminosity (e.g., Shen et al. 2008a; Fine et al. 2008), i.e., for any given luminosity, the FWHM distribution has a constant mean value and a log-normal scatter σ_FWHM around the mean. Therefore at fixed luminosity, the value of m_e is independent of the true mass m, i.e.,

$$\begin{equation} f(m_e|\lambda)=\big(2\pi \sigma _{\rm line}^2\big)^{-1/2}\exp \left\lbrace -\frac{[m_e-(C_1\lambda +C_2)]^2}{2\sigma _{\rm line}^2}\right\rbrace {,} \end{equation} \tag{ 10 }$$

where σ_line = 2σ_FWHM is the scatter in virial mass resulting from the scatter in FWHM (σ_FWHM) at fixed luminosity, and the constant mean value of FWHM is absorbed in C₂.

Equations (6)–(8) have analytical results for specific forms of Ψ_M(m) and g(λ|m). For a power-law true BH mass function $\Psi _M(m)\propto 10^{m\gamma _M}$, and a power-law relation between true mass and luminosity 〈λ〉 = am + b with a constant log-normal scatter, σ_L, i.e.,

$$\begin{equation} g(\lambda |m)=\big(2\pi \sigma _L^2\big)^{-1/2}\exp \left\lbrace -\frac{[\lambda -(am+b)]^2}{2\sigma _L^2}\right\rbrace {,} \end{equation} \tag{ 11 }$$

Equations (6)–(8) yield

$$\begin{equation} \Psi _L(\lambda)\propto 10^{\lambda \gamma _M/a} {;} \end{equation} \tag{ 12 }$$

$$\begin{eqnarray} \langle m\rangle &=&\frac{\gamma _M\ln (10)\sigma _L^2-ab}{a^2}-\frac{1}{\gamma _M\ln (10)}\nonumber \\ &&+\frac{\lambda _2 10^{\lambda _2\gamma _M/a}-\lambda _1 10^{\lambda _1\gamma _M/a}}{a(10^{\lambda _2\gamma _M/a}-10^{\lambda _1\gamma _M/a})}{;} \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \langle m_e\rangle &= &C_2 -\frac{C_1a}{\gamma _M\ln (10)}\nonumber \\ &&+\frac{C_1(\lambda _2 10^{\lambda _2\gamma _M/a}-\lambda _1 10^{\lambda _1\gamma _M/a})}{10^{\lambda _2\gamma _M/a}-10^{\lambda _1\gamma _M/a}} {.} \end{eqnarray} \tag{ 14 }$$

Note throughout this paper we assume that this constant scatter in luminosity, σ_L, is solely the variations that are uncompensated by the corresponding variations in FWHM, i.e., it is one of the two sources of the virial uncertainty σ_vir (Equation (1)).³

There are eight parameters that determine the bias Δlog M_• = 〈m_e〉 − 〈m〉: σ_L, γ_M, a and b are for the underlying true BH mass function and luminosity (or Eddington ratio) distributions at fixed true mass, which are model assumptions and must be tuned to produce the observed quasar LF; C₁ and C₂ are determined from the empirically calibrated virial estimators and measured FWHM distributions in statistical quasar samples; λ₁ and λ₂ are observational windows determined from the specific quasar sample under investigation. A further constraint comes from the assumption that virial mass is an unbiased estimator of the true BH mass (e.g., Equation (2)), which requires (e.g., consulting Equations (2), (10) and (11))

$$\begin{equation} C_1=1/a,\qquad C_2=-b/a{.} \end{equation} \tag{ 15 }$$

Thus the bias reduces to

$$\begin{equation} \Delta \log M_{\bullet }=-\gamma _M\ln (10)\sigma _L^2\big/a^2{,} \end{equation} \tag{ 16 }$$

from Equations (13) and (14) (see also Section 4.4 in Shen et al. 2008a), which is independent of luminosity.

For demonstration purposes here we take a model with parameters σ_L = 0.4, γ_M = −4, a = 2, and b = 28.4. These parameters produce a slope of γ_L = −2 in the quasar LF, close to the observed bright-end slope (e.g., Richards et al. 2006b; Hopkins et al. 2007), and constants C₁ = 1/a = 0.5 and C₂ = −b/a = −14.2, which are consistent with the commonly used virial mass calibrations and the observed FWHM distributions (e.g., McLure & Jarvis 2002; Vestergaard & Peterson 2006; McGill et al. 2008; Shen et al. 2008a). With these parameters, the difference between the sample-averaged true and virial BH masses is 〈m_e〉 − 〈m〉 ≈ 0.37 dex.

The derivation of this virial mass bias is independent from the derivation of the bias discussed in Lauer et al. (2007). To illustrate this, let us consider again simple power-law distributions and Gaussian (log-normal) scatter. We start from a power-law distribution of quasar host galaxy property s (s ≡ log M_bulge for instance),

$$\begin{equation} \Psi _S(s)\propto 10^{s\gamma _S}{.} \end{equation} \tag{ 17 }$$

The distribution of true BH mass m at fixed s is

$$\begin{equation} p_0(m|s)=\big(2\pi \sigma _{\mu }^2\big)^{-1/2}\exp \left\lbrace -\frac{[m-(C_3s+C_4)]^2}{2\sigma _{\mu }^2}\right\rbrace {,} \end{equation} \tag{ 18 }$$

where the BH scaling relation is 〈m〉 = C₃s + C₄ with a cosmic scatter σ_μ, where C₃ and C₄ are constants. Given the power-law distribution of s (Equation (17)), the distribution of true BH mass m is then also a power-law $\Psi _M(m)\propto 10^{m\gamma _M}$ where γ_M = γ_S/C₃. Assuming the distribution of luminosity λ at fixed m is given by Equation (11) we can derive the distribution of λ at fixed s as

$$\begin{eqnarray} &&p_0(\lambda |s)=\int g(\lambda |m)p_0(m|s) dm\nonumber \\ &&\quad =\big[2\pi \big(\sigma _L^2+a^2\sigma _{\mu }^2\big)\big]^{-1/2}\exp \left\lbrace -\frac{\displaystyle [\lambda -(aC_3s+b+aC_4)]^2}{\displaystyle 2\big(\sigma _L^2+a^2\sigma _{\mu }^2\big)}\right\rbrace. \nonumber \\ \end{eqnarray} \tag{ 19 }$$

Therefore the distribution of s at fixed luminosity λ is

$$\begin{eqnarray} p_1(s|\lambda)&=&p_0(\lambda |s)\Psi _S(s)\left(\int p_0(\lambda |s)\Psi _S(s)ds\right)^{-1}\nonumber \\ &=&\big(2\pi \sigma _1^2\big)^{-1/2}\exp \left\lbrace -\frac{[s-(C_5\lambda +C_6)]^2}{2\sigma _1^2}\right\rbrace, \end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray} &&\sigma _1^2=\frac{\sigma _L^2+a^2\sigma _\mu ^2}{a^2C_3^2},\ C_5=\frac{1}{aC_3},\ C_6=\ln (10)\gamma _S\sigma _1^2 - \frac{b+aC_4}{aC_3}{.}\nonumber \\ \end{eqnarray} \tag{ 21 }$$

Similarly, the distribution of m at fixed λ is

$$\begin{eqnarray} p_1(m|\lambda)&=&g(\lambda |m)\Psi _M(m)\left(\int g(\lambda |m)\Psi _M(m)dm\right)^{-1}\nonumber \\ &=&\big(2\pi \sigma _2^2\big)^{-1/2}\exp \left\lbrace -\frac{[m-(C_7\lambda +C_8)]^2}{2\sigma _2^2}\right\rbrace\! {,} \end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray} &&\sigma _2^2=\frac{\sigma _L^2}{a^2},\ C_7=\frac{1}{a},\ C_8=\frac{\ln (10)\gamma _S\sigma _L^2}{a^2C_3}-\frac{b}{a}{.} \end{eqnarray} \tag{ 23 }$$

Hence, the Lauer et al. bias, i.e., the excess in true BH mass at fixed luminosity is simply

$$\begin{eqnarray} \Delta \log M_{\bullet,{\rm Lauer}}&=&\langle m \rangle |_\lambda - (C_3\langle s \rangle |_\lambda + C_4)\nonumber \\ &=&{-}\ln (10)\gamma _S\sigma _\mu ^2\big/C_3{.} \end{eqnarray} \tag{ 24 }$$

Since γ_S is negative, the Lauer et al. bias states that at fixed luminosity, the average true BH mass is already biased high with respect to the expectation from the mean BH scaling relation. On top of that, the virial mass bias (Equation (16)) states that the average BH mass estimate is further biased high with respect to the true mass.

The simple estimate of the BH mass bias discussed above neglects the curvature in the intrinsic BHMF and quasar LF, and hence is only valid at the bright end of the LF. In particular, since the quasar LF flattens below the break luminosity λ_break ∼ 46, we expect that the mass bias becomes less severe toward fainter luminosities. It is important to realize that in essentially all of the studies mentioned in Section 1, larger offsets usually occur at higher redshift, when their samples are sampling the higher-mass end of the intrinsic BHMF. It is conceivable then that a false evolution will arise simply from the increasing mass bias toward the high-mass tail.

To estimate the BH mass bias at various redshifts and luminosity sampling ranges, we assume an underlying BHMF and a model for the luminosity distribution at fixed true BH mass,⁴ using the observed quasar LF as constraints. This forward-modeling procedure is similar to the simulations done in Shen et al. (2008a) and Kelly et al. (2009a, 2009b).

We start with a broken power-law model for the true BHMF:

$$\begin{equation} \Psi _{M}(m)\propto \frac{1}{10^{-\gamma _{M1}(m-m_*)}+10^{-\gamma _{M2}(m-m_*)}}\;, \end{equation} \tag{ 25 }$$

where γ_M1 and γ_M2 are the low-mass and high-mass end slope, and m_* is the break BH mass. For the luminosity distribution g(λ|m), we still use the Gaussian form in Equation (11). For simplicity, we assume that the set of six parameters, γ_M1, γ_M2, m_*, a, b, and σ_L, are only functions of redshift. At a fixed redshift, a set of the six parameters determines the shape of the quasar LF (7), which is to be constrained by observations. Assuming that the commonly used virial BH estimators are unbiased, we fix a = 2 and b = 28.4, which corresponds to a BLR radius–luminosity relation R ∝ L^0.5 and virial coefficients and FWHM distributions consistent with previous work (e.g., Shen et al. 2008a; Fine et al. 2008). We also fix the value of σ_L during model fitting. Therefore, there are four free parameters: γ_M1, γ_M2, m_*, and the normalization of the BHMF, Φ₀ (in units of Mpc⁻³dex⁻¹). We list the best-fit model parameters in Table 1 for different values of σ_L and at several redshifts. We note that these models are mainly used to demonstrate the effects of the virial mass bias and should not be interpreted as our attempt to constrain the true active BHMF.

Figure 1 shows the model LF with fixed σ_L = 0.4 and comparison with the bolometric LF data compiled in Hopkins et al. (2007, and references therein) at several redshifts. At each redshift, we determine γ_M1, γ_M2, m_*, and Φ₀ by minimizing the χ² between the model and the data. As discussed earlier, this scatter σ_L refers to the variations that are uncompensated by the corresponding variations in FWHM. If we were to incorporate an additional correlated variation term in luminosity at fixed true mass, which are compensated by an additional variation in FWHM, we would need a steeper high-mass end slope for the true BHMF in order not to overshoot the bright-end LF, which then makes the mass bias even worse. Also, the amount of this additional broadening is limited by the already narrow distributions of FWHM (e.g., Shen et al. 2008a; Fine et al. 2008) and the requirement that σ_vir ≳ 0.4.

Zoom In Zoom Out Reset image size

Given the model true BHMF and luminosity distribution, we can now estimate the bias Δlog M_• = 〈m_e〉 − 〈m〉 as a function of the sampled luminosity range using Equations (6) and (8). Figure 2 shows Δlog M_• for luminosity-limited quasar samples with λ>λ₀ ≡ log L₀ at several redshifts and for fixed σ_L = 0.3, 0.4, and 0.5 dex. The bias generally rises when luminosity increases and approaches the asymptotic value −γ_M2ln(10)σ²_L/a² at the bright end. Because we are fitting our model against the bolometric LF data at individual redshifts, we do not expect identical results since there will be systematics involved in deriving the bolometric LF data from band conversions and from the assumed quasar spectral energy distribution (see details in Hopkins et al. 2007). Nevertheless, this analysis demonstrates that the mass bias 〈m_e〉 − 〈m〉 could be substantial for bright luminosities and for large dispersion σ_L in the luminosity distribution at fixed true mass.

Zoom In Zoom Out Reset image size

The virial mass bias discussed here depends on the amplitude of the scatter σ_L. This is the variation in luminosity at fixed true BH mass which is not compensated by the variation in FWHM. How large is σ_L? For the best-studied local RM samples and for Hβ only (e.g., Kaspi et al. 2005; Bentz et al. 2006, 2009), the scatter of luminosity around fixed BLR size is on the order of ≲0.2 dex. However, this is for the best cases. We generally expect larger σ_L for statistical quasar samples with single-epoch spectra for the following reasons: (1) single-epoch spectra have larger scatter than averaged spectra in RM samples; (2) these samples usually cover a luminosity range that is poorly sampled by current RM data; (3) the UV regime of the spectrum (in particular for C iv) has more peculiarities and larger scatter in the R–L relation than the Hβ regime, but is usually used at high redshift due to the drop off of Hβ in the optical spectrum; (4) even if some variations in luminosity truthfully trace the variations in BLR radius, R, they may still be uncompensated by FWHM if, for instance, there is a non-virial component in the broad line which does not respond to BLR radius variations (especially for C iv; see discussions in Shen et al. 2008a); (5) in many statistical quasar samples, the spectra have low signal-to-noise ratios, which lead to increased uncorrelated scatter between the measured luminosity and FWHM, and bias the virial mass estimates more. Hence it is very plausible that σ_L ≳ 0.4 in most cases, which then accounts for ≳0.2 dex in σ_vir (e.g., Equation (1)).

Therefore, we conclude that the mass bias discussed here contributes at least 0.2–0.3 dex at L_bol ≳ 10⁴⁶ erg s^-1, comparable to the statistical bias resulting from the cosmic scatter in the BH scaling relations (e.g., Lauer et al. 2007). In Figure 2, we also show the "observed" BH mass offset (excess at fixed galaxy properties) from the literature. We note that the amount of the uncorrelated scatter σ_L might increase with redshift both due to the switch from the Hβ line to the more problematic UV lines and due to the often decreased spectral quality at high redshift. Since the virial mass bias is independent of the Lauer et al. bias, the two statistical biases together can account for a large (or even full) portion of the BH mass excess seen in the data, hence no exotic scenarios are needed to explain this strong redshift evolution. It is possible, however, that there is still a mild evolution in these BH scaling relations, which is inherent to the cosmic co-evolution of SMBHs and their hosts. But given these statistical biases, and given other systematics with single-epoch virial estimators (not discussed here, but see, e.g., Denney et al. 2009), it is premature to claim a strong positive evolution in the BH scaling relations.

We thank the anonymous referee, as well as Tod Lauer, Scott Tremaine, Jenny Greene, and Avi Loeb for useful comments. Y.S. acknowledges support from a Clay Postdoctoral Fellowship through the Smithsonian Astrophysical Observatory (SAO). B.K. acknowledges support by NASA through Hubble Fellowship grant No. HF-51243.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.