The Mass Relations between Supermassive Black Holes and Their Host Galaxies at 1 < z < 2 with HST-WFC3

The 32 AGNs were initially selected by the X-ray observations of COSMOS (Civano et al. 2016), (E)-CDFS-S (Lehmer et al. 2005; Xue et al. 2011), and SXDS (Ueda et al. 2008) fields. In most cases, the X-ray sources are first identified as broad-line AGNs through optical spectroscopic campaigns with the Very Large Telescope, Keck, and Magellan. Follow-up near-infrared spectroscopic observations of the AGNs in these fields are carried out with Subaru's FMOS (Kimura et al. 2010; Nobuta et al. 2012; Matsuoka et al. 2013), covering the wavelength range 0.9−1.8 μm, which provides the favorable Hα and Hβ lines out to z ∼ 1.7 to estimate the ${{ \mathcal M }}_{\mathrm{BH}}$. Recently, Schulze et al. (2018) presented near-IR spectroscopy of a large compilation of 243 X-ray AGNs in these fields. It is from this catalog that we primarily select our targets. The continuum fitting and emission-line modeling provided by this work are performed through a common procedure for spectral model fits. This approach first corrects the spectra for galactic extinction and shifts them to the rest frame. The spectral regions in the near-IR that are strongly affected by OH emission are masked. The spectrum is modeled by both a narrow-line and a broad-line AGN template, and the emission-line information is inferred based on χ² minimization. We refer the reader to the aforementioned paper for further details. Note that three CDFS objects (i.e., CDFS-1, CDFS-229, and CDFS-724) are not included in Schulze et al. (2018); their infrared spectra are provided by Suh et al. (2015) using a similar approach.

Based on the ${{ \mathcal M }}_{\mathrm{BH}}$ estimates described below (Section 2.2), we select targets with masses in the range 7.5 ≲ log (${{ \mathcal M }}_{\mathrm{BH}}$/M_⊙) ≲ 8.5. The bolometric luminosities are measured from the broad AGN lines by Schulze et al. (2018; Section 3.3), which are used to calculate their Eddington ratio λ = L_bol/L_Edd. The flux-limited nature of the sample results in a slightly higher Eddington ratio distribution at lower ${{ \mathcal M }}_{\mathrm{BH}}$, as shown in Figure 1. Moreover, we have a higher preference to select targets that have rest-frame UV images (i.e., HST/ACS), as provided by Scoville et al. (2007) and Koekemoer et al. (2007) in the COSMOS field. We list the 32 AGNs observed with HST/WFC3 and analyzed in this work in Table 1, sorted by field and redshift.

The ${{ \mathcal M }}_{\mathrm{BH}}$ of type 1 AGNs can be determined using the so-called virial method (Peterson et al. 2004; Shen 2013). The kinematics of the BLR trace the gravitational field of the central SMBH, assuming the gravity dominates the motion of the BLR gas. Under these assumptions, the width of the emission line provides the scale of the velocity dispersion (ΔV), while the AGN continuum luminosity establishes an empirical scale of the BLR size (R_BLR). The estimation of ${{ \mathcal M }}_{\mathrm{BH}}$ is then achieved using these measurements, i.e., ${{ \mathcal M }}_{\mathrm{BH}}$ ≃ G⁻¹ R_BLR ΔV² (McLure & Dunlop 2004).

To avoid any systematic bias between samples in the literature, we adopt a class of self-consistent estimators for our analysis. We first compare the estimators implemented in Schulze et al. (2018) and Ding et al. (2017b) and find very consistent Hβ(FWHM(5100))-based masses (${{ \mathcal M }}_{\mathrm{BH}}$ rms <0.03 dex). However, there is an ∼0.2 dex inconsistency in their Hα mass estimates. Therefore, we utilize the AGNs in Schulze et al. (2018) that have both Hα and Hβ lines (35 AGNs in total) and carry out a cross-calibration to determine which Hα estimator has the best agreement between the two lines. As a result, the Hα estimator in Schulze et al. (2018) has better agreement with Hβ. Thus, we adopt the scheme given in Schulze et al. (2018) for all AGN samples used in this study, including the comparison samples described below:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\left(\displaystyle \frac{{{ \mathcal M }}_{\mathrm{BH}}({\rm{H}}\alpha )}{{M}_{\odot }}\right) & = & 6.71+0.48\mathrm{log}\left(\displaystyle \frac{{L}_{{\rm{H}}\alpha }}{{10}^{42}\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)\\ & & +\,2.12\mathrm{log}\left(\displaystyle \frac{\mathrm{FWHM}({\rm{H}}\alpha )}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\end{array}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\left(\displaystyle \frac{{{ \mathcal M }}_{\mathrm{BH}}({\rm{H}}\beta )}{{M}_{\odot }}\right) & = & 6.91+0.50\mathrm{log}\left(\displaystyle \frac{{L}_{{\lambda }_{5100}}}{{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)\\ & & +\,2.0\mathrm{log}\left(\displaystyle \frac{\mathrm{FWHM}({\rm{H}}\beta )}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right).\end{array}\end{eqnarray} \tag{ 2 }$$

Note that these recipes are first provided by Vestergaard & Peterson (2006). Using these recipes, we estimate ${{ \mathcal M }}_{\mathrm{BH}}$ by adopting the emission-line properties as measured by Schulze et al. (2018) for all 32 AGNs. While 14 AGNs have emission-line properties for both Hα and Hβ, we adopt the value of ${{ \mathcal M }}_{\mathrm{BH}}$ based on the Hα emission line to have consistency across the sample. While most recipes are calibrated against Hβ, the line is typically weaker than Hα, hence lower signal-to-noise. We provide the ${{ \mathcal M }}_{\mathrm{BH}}$ measurements, together with the properties of the emission lines, in Table 2.

It is worth noting that the absolute flux calibration of the FMOS spectra is set to match the available ground-based IR imaging, UltraVISTA in the case of COSMOS. While an initial flux calibration is performed during the reduction of the FMOS data using calibration stars, there can be differential flux loss due to aperture effects, variable seeing conditions, and minor alignment issues with the instrument. The flux normalization is effectively an aperture correction in the J or H band, depending on the source redshift. We note that this procedure does not correct for intrinsic variability that may induce an additional error of 0.2 mag. We refer the reader to Section 2.2 of Schulze et al. (2018) for full details of the flux calibration.

Our AGN sample is primarily selected based on the value of ${{ \mathcal M }}_{\mathrm{BH}}$ and the Eddington ratio. It is well known that sample selection effects must be taken into account in order to interpret the observed black hole–host relations and measure their evolution with redshift, thus avoiding biases (Treu et al. 2007; Bennert et al. 2011a; Schulze & Wisotzki 2011, 2014; Park et al. 2015). The main source of bias is due to the fact that active samples are necessarily selected based on properties that correlate with black hole mass estimators (such as AGN luminosity, line strength, and width), and thus one tends to favor overly massive black holes in the presence of intrinsic scatter and observational errors. Correcting for observational biases requires a well-characterized selection function, such as the one we have for our sample; a model of the black hole mass function; and the evolution of the correlations between ${{ \mathcal M }}_{\mathrm{BH}}$ and other properties.

Given our selection function, we use the Bayesian framework introduced by Schulze & Wisotzki (2011, 2014) to estimate the expected bias. In this framework, under the assumption of no evolution of the correlations between ${{ \mathcal M }}_{\mathrm{BH}}$ and host galaxy properties, one can compute the expected bias for a given sample prior to the observations. The key ingredients of this model are the local ${{ \mathcal M }}_{\mathrm{BH}}$ host galaxy property correlations, black hole mass function, and Eddington ratio distribution at the redshift of observation. The latter two quantities are estimated from the type 1 AGN distributions (Schulze et al. 2015) and corrected to represent the parent population of all active SMBH hosts.

Adopting our specific selection limits into the framework (i.e., $\mathrm{log}($ ${{ \mathcal M }}_{\mathrm{BH}}$ $/{M}_{\odot })\in [7.5,8.56]$, $\mathrm{log}({L}_{\mathrm{bol}}/{L}_{\odot })\in [45.0,46.2]$, and $\mathrm{log}(\lambda )\in [-2.0,0.5]$), we infer an expected bias of +0.21 dex in the observed ${\rm{\Delta }}\mathrm{log}({{ \mathcal M }}_{\mathrm{BH}}$) for samples at z ∼ 1.5, assuming the baseline choices for the inputs to the model that include the local value of the mass ${{ \mathcal M }}_{\mathrm{BH}}$/M_* ratio. To reiterate, an offset in the observed mass relations of 0.2 dex at z ∼ 1.5 can be considered consistent with no evolution in the mass ratio. This bias correction should be considered an approximate estimate, with an uncertainty that depends on the uncertainty of the inputs. Furthermore, if the scaling relations actually evolve, one needs to introduce a model for the evolution in order to infer the underlying bias-corrected trends. We will revisit these issues in Section 5.5.

We make use of the ${{ \mathcal M }}_{\mathrm{BH}}$–L_host and ${{ \mathcal M }}_{\mathrm{BH}}$–M_* relations in the literature for comparison with our high-z data. For the local relation, we use the measurements of Bennert et al. (2010) and Bennert et al. (2011a, hereafter B10 and B11) to define our zero-point for local AGN samples. The sample by B10 consists of 19 AGNs with the ${{ \mathcal M }}_{\mathrm{BH}}$–L_host relation determined with reliable ${{ \mathcal M }}_{\mathrm{BH}}$ masses using reverberation mapping with an uncertainty level of ∼0.15 dex. Note that B10 only provided the single-galaxy V-band luminosity. Ding et al. (2017b) derived the galaxy R-band luminosity based on the same early-type galaxy template spectrum adopted by B10 using K-correction (see Section 3.2 therein). The work of B11 contains 25 local active AGNs, where the ${{ \mathcal M }}_{\mathrm{BH}}$ are measured using the single-epoch method (${{ \mathcal M }}_{\mathrm{BH}}$ uncertainty level ∼0.4 dex). To track the local M_* and ${{ \mathcal M }}_{\mathrm{BH}}$ relations to higher values, we include 30 inactive galaxies (mainly ellipticals or S0) from Häring & Rix (2004, hereafter HR04). It is worth noting that the local inactive sample is mainly bulge-dominated, and we adopt the bulge mass for the entire local sample. In other words, our local comparison is the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,bulge and not those involving total quantities.

We include in our analysis published samples at intermediate redshifts to understand the evolution of these correlations. We select samples that were analyzed by members of our team to ensure uniform measurements. The intermediate-redshift AGNs that we include are 52 objects published by Park et al. (2015) using a single band that are applicable for the ${{ \mathcal M }}_{\mathrm{BH}}$–L_host relation at 0.36 < z < 0.57 and 27 objects published by Bennert et al. (2011b) and Schramm & Silverman (2013) at 0.5 < z < 1.9. Similar to the B10 sample, the R-band luminosities of these intermediate-redshift systems (79 in total) are obtained by Ding et al. (2017b) using K-corrections from the V band based on the same stellar template. Moreover, Bennert et al. (2011b) and Schramm & Silverman (2013) estimated the stellar masses of their 27 objects using multiband imaging data, and we adopt them as the ${{ \mathcal M }}_{\mathrm{BH}}$–M_* comparison sample. In addition, we adopt a sample of 32 ${{ \mathcal M }}_{\mathrm{BH}}$–M_* measurements at 0.3 < z < 0.9 by Cisternas et al. (2011) to compare with our sample. Across the intermediate-redshift sample, we recalibrate the ${{ \mathcal M }}_{\mathrm{BH}}$ using the self-consistent recipes introduced in Section 2.2 including those estimated from Mg_II using the recipe in Ding et al. (2017b).

The values for all of the ${{ \mathcal M }}_{\mathrm{BH}}$–L_host comparison samples are listed in Ding et al. (2017b; Tables 1 and 2 therein). The ${{ \mathcal M }}_{\mathrm{BH}}$–M_* comparison samples are collected and summarized here in Appendix B, Table 7.

The knowledge of the PSF is crucial for imaging the decomposition of the AGN and its host, especially when the point source contributes to the majority of the total emission. The PSF is known to vary across the detector and over time due to the effects of aberration and breathing. Simulated PSFs, such as those based on TinyTim, are usually insufficient for our purposes (Mechtley et al. 2012). Stars within the field of view of each observation provide a better description than the simulated PSF, since they are observed simultaneously with the science targets and reduced and analyzed in a consistent manner (Kim et al. 2008; Park et al. 2015). However, we have found that such stars usually do not provide an ideal PSF for deblending the AGN and host galaxy through extensive tests. The issues are associated with an insufficient number of bright stars near our targets, color differences, and other effects not fully understood.

To minimize the impact of such mismatches, we build a PSF library by selecting all of the isolated, unsaturated PSF stars with high signal-to-noise ratio from our entire program. The selection consists of the following steps. First, we identify stars from the COSMOS2015 catalog (Laigle et al. 2016). However, many bright stars with an intensity similar to our AGN sample were excluded in this catalog. Therefore, we also manually select PSF-like objects as candidates from the HST/WFC3 imaging. We then discard nonideal PSF candidates based on their intensity, FWHM, central symmetry, and presence of nearby contaminants. In total, the PSF library contains 78 and 37 stars imaged through filters F140W and F125W, respectively. We assume that the stars in the library are representative of the possible PSFs in our program. The dispersion of PSF shapes within our library provides us with a good representation of the level of uncertainty in our measurements resulting from the image decomposition.

As shown in Figure 3 (top panel), the host-to-total flux ratio (total = host+nuclear) of the sample spans a wide range, from 10% to 90%, with most of the sample concentrated between 20% and 50% (median value 37%). Those with flux ratios above ∼20% have a higher degree of significance with respect to the detection of the host galaxy, while the five systems (i.e., CID 255, CID 50, LID 360, CDFS-229, SXDS-X763) that have host-to-total flux ratios lower than 20% should be considered as marginal detections.

The distributions of effective radius R_eff are shown in Figure 3 (bottom panel). Based on the redshift, we calculate the physical scale of the radius for each object in kpc assuming a standard cosmology. We plot them together with the measured Sérsic index in Figure 4. These values are distributed within the allowed range and not concentrated on either the upper or lower bound. The R_eff values of our sample are between ∼1 and 7 kpc (peaked at ∼2.2 kpc). Nearly half (15/32) of our systems have a Sérsic index n < 2, indicating that they have a significant disk component, likely in addition to the presence of a bulge. Five systems have a high Sérsic index (n > 4.5) for the first run; we use n ∈ [1, 4] as a prior to refit these systems and find that the changes on the inference of their host luminosity are very limited (<0.03 dex). Seven systems have an effective radius R_eff < 018, i.e., smaller than 3 pixels. In order to check whether this affects our conclusions, we refit them with an R_eff ∈[02, 10] prior and find that the inferred host flux barely changes (<1%). In particular, the inferred R_eff for three systems (CID 543, XID 2202, and SXDS-X1136) hit the lower limit (i.e., 01); thus, the scatter on the parameter is formally zero. This means that the size inference of these three systems should really be considered an upper limit, rather than a measurement. However, the inferences of their host luminosities and stellar masses are still reliable given the consistency of the refitting results and the small scatter in the host flux (<10%).

Zoom In Zoom Out Reset image size

We compare the morphology of AGN host galaxies to inactive galaxies from the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). We identify 4401 inactive galaxies within a redshift range (1.2 < z < 1.7), comparable to our AGN sample, whose Sérsic measurements are provided by van der Wel et al. (2012) using Galfit. Their stellar masses are derived based on the 3D-HST spectroscopic survey (Brammer et al. 2012; Momcheva et al. 2016) and comparable to our AGN hosts ($9.5\lt \mathrm{log}({M}_{* }/{M}_{\odot })\lt 11.5;$ Section 5.4). We compare the histogram of the inferred R_eff and Sérsic index to the inactive galaxies in Figure 4, where we find no significant difference between their distributions and median value. We also test whether the distributions can be drawn from the same parent population using a Kolmogorov–Smirnov test and determining the p-value to be 0.42 and 0.04 for R_eff and n, respectively. We conclude that the host galaxies of our AGN sample are representative of the overall population of galaxies, with a significant disk component, at comparable luminosity and stellar mass at the same redshift. In Section 5.6, we use this information to infer the likely bulge masses, hence the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,bulge relation.

For 21/32 AGNs, we have multiband host magnitudes that enable us to select the appropriate stellar population templates to determine rest-frame luminosities and stellar masses for the overall sample. We find that the 1  and 0.625 Gyr stellar populations with solar metallicity and a Chabrier IMF (Bruzual & Charlot 2003) provide excellent matches to the observed colors of our sample at z < 1.44 and z > 1.44, respectively (see Figure 5). To minimize the uncertainty associated with these corrections, we use these two templates to interpolate to the rest-frame R band, which is very close to the observed wavelengths. We note that the choice is not unique, and other combinations of ages, metallicities, and star formation history could match the observed colors and would provide very similar R-band magnitudes, as well as stellar masses (Bell & de Jong 2000, 2001).

Zoom In Zoom Out Reset image size

Adopting the 1  and 0.625 Gyr stellar populations, we perform a K-correction to derive the rest-frame R-band magnitude of our sample, based on the host inference in the WFC3 band. As mentioned above, since the WFC3 filter is already close to the rest-frame R band, we expect the M_R uncertainty introduced by this K-correction to be within 0.05 mag. We derive the rest-frame R-band luminosity from $\mathrm{log}{L}_{R}/{L}_{R,\odot }=0.4\times ({M}_{R,\odot }\mbox{--}{M}_{R})$, where ${M}_{R,\odot }=4.61$ (Blanton & Roweis 2007). Here L_R ranges between $\mathrm{log}({L}_{R}/{L}_{\odot ,R})\in [9.5,11.5]$ with individual values listed in Table 4.

We show the relation between ${{ \mathcal M }}_{\mathrm{BH}}$ and L_host in Figure 6 with comparison samples. For reference, we fit the local data with a linear relation,

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{{ \mathcal M }}_{\mathrm{BH}}}{{10}^{7}{M}_{\odot }}\right)={\alpha }_{0}+{\beta }_{0}\mathrm{log}\left(\displaystyle \frac{{L}_{R}}{{10}^{10}{L}_{\odot }}\right),\end{eqnarray} \tag{ 7 }$$

which enables a direct comparison between our high-z sample and the local relation. The distribution of our data appears to be in good agreement with the local relation and the other AGN samples at lower redshift. Therefore, the observational data indicate that the relation between black hole mass and host luminosity are similar at different periods of the universe. In Appendix A, we explore how the high-z sample would evolve in this plane solely with the luminosity evolution of the host galaxy as done in past studies (e.g., Ding et al. 2017b).

Zoom In Zoom Out Reset image size

We note that there are outliers that deviate from the distribution of the overall sample, such as SXDS-X763. We suspect that this source may have abnormal host properties. Indeed, the value of R_eff has a large uncertainty (∼75%), and the host-to-total flux ratio is the lowest of the high-z sample (<10%).

Using the near-IR imaging with HST, we take the host luminosity along with color information to estimate the stellar mass content of each host galaxy based on the mass-to-light ratio of the adopted stellar populations. The uncertainty level associated with the stellar mass is expected to be of order 0.1 dex (changing the IMF would affect all of our stellar masses systematically). We find that M_* ranges between $\mathrm{log}({M}_{* }/{M}_{\odot })\in [9.7,11.3]$. These values are listed in Table 4.

In Figure 7, we plot our ${{ \mathcal M }}_{\mathrm{BH}}$–M_* measurements and find that the distributions in this plane for the four samples are similar. However, there is a slight shift of the high-z data toward higher black hole masses at a fixed host mass. In Figure 8, we plot the mass ratio as a function of redshift (right panel) and the mass ratio difference (${\rm{\Delta }}\mathrm{log}{{ \mathcal M }}_{\mathrm{BH}}=\mathrm{log}\left(\tfrac{{{ \mathcal M }}_{\mathrm{BH}}}{{10}^{7}{M}_{\odot }}\right)-\alpha -\beta \mathrm{log}\left(\tfrac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)$, which is considered as the difference between the high-redshift data and the best-fit local relation; left panel). For our high-z sample, we find the averaged ${\rm{\Delta }}\mathrm{log}{{ \mathcal M }}_{\mathrm{BH}}$ to be 0.43 ± 0.06, a factor of ∼2.7 higher than the local mass ratio, as indicated by the top blue circle in Figure 8.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As with Equation (7), we fit the ${{ \mathcal M }}_{\mathrm{BH}}$–M_* data with a linear relation as

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{{ \mathcal M }}_{\mathrm{BH}}}{{10}^{7}{M}_{\odot }}\right)={\alpha }_{1}+{\beta }_{1}\mathrm{log}\left(\displaystyle \frac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)\end{eqnarray} \tag{ 8 }$$

and an evolution term parameterized in the following form:

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{log}{{ \mathcal M }}_{\mathrm{BH}}=\gamma \mathrm{log}(1+z).\end{eqnarray} \tag{ 9 }$$

Performing such a fit based on the observed data (i.e., no correction for selection effects), we obtain γ = 1.03 ± 0.25. Even with the exclusion of the few outliers, the slope is significantly nonzero with γ = 0.89 ± 0.27. The best-fit observed evolution model is shown in Figure 8.

As shown in these panels, the scatter of the high-z correlation presents a similarity to the local relation. To quantitatively make this comparison, we investigate the intrinsic scatter of our sample. Note that the observed high-z sample has effects from both the measurement uncertainty and the narrow selection window, and thus one needs to take them both into account to extract the intrinsic scatter. Based on our forward-modeling framework, we find that our AGN sample has intrinsic scatter as 0.25 dex on the vertical axis, taking into account the measurement uncertainty and the selection function. We will return to this topic in a forthcoming paper (X. Ding et al. 2019, in preparation), where we focus on the scatter of the sample as a diagnostic of the AGN–host galaxy feedback mechanism. The fact that the intrinsic scatter for the high-z sample is no larger than the local one, which is ∼0.35 dex (Gültekin et al. 2009), may pose a challenge for explanations of scaling relations where random mergers are the origin of the correlation (Peng 2007; Jahnke & Macciò 2011), indicating that there may likely be a connection between the SMBHs and their host galaxies during their formation.

In Section 2.3, we used the Bayesian framework introduced by Schulze & Wisotzki (2011) and found that in the case of no intrinsic evolution of the correlations, we would expect to measure an offset corresponding to a +0.21 dex bias in the inferred ${\rm{\Delta }}\mathrm{log}$ ${{ \mathcal M }}_{\mathrm{BH}}$, owing to our selection function. Interestingly, this bias could account for the majority of the observed offset as obtained in the ${{ \mathcal M }}_{\mathrm{BH}}$–M_* correlations in the last section. To demonstrate, we use blue open circles to show the expected bias raised by this selection effect in Figure 8. We also show, in Figure 9, the histogram ${\rm{\Delta }}\mathrm{log}$ ${{ \mathcal M }}_{\mathrm{BH}}$ of the high-z sample and compare it to the local ones. The plots show that the correlations between ${{ \mathcal M }}_{\mathrm{BH}}$ and host galaxy total stellar mass or luminosity could be consistent with those of the local samples, given the uncertainties and selection function.

Zoom In Zoom Out Reset image size

To further evaluate the "true" underlying evolution and its uncertainty, corrected for selection effects and taking into account the actual observations, we adopt an independent method based on the approach introduced by Treu et al. (2007) and developed by Bennert et al. (2010), Park et al. (2015), and Ding et al. (2017b). The method parameterizes the evolution of the correlations between ${{ \mathcal M }}_{\mathrm{BH}}$ and host properties as an offset from the local γ and an intrinsic scatter σ_int, which can be a free parameter or tied to the local relation. It then imposes a selection function in ${{ \mathcal M }}_{\mathrm{BH}}$ and calculates the intrinsic parameters of the model given the observations. In practice, we start from the local black hole mass function and the evolution model, generate mock samples using a Monte Carlo approach, apply the selection function, and compare with the data to generate the likelihood. To illustrate the importance of the intrinsic scatter in the selection effects, we adopt both a uniform (flat) prior and a lognormal prior for σ_int. Note that this method assumes a narrow Eddington ratio distribution, which is different from the one by Schulze & Wisotzki (2011), as described in Section 2.3. Also, this method adopts the local black hole mass function rather than the high-z black hole mass function. These differences have a second-order effect and could be responsible for the different magnitude of the selection effect. On the other hand, this method probes the importance of the scatter at high-z, which is complementary.

Combining the 32 AGNs together with the intermediate-redshift sample, we present the inferred γ and σ_int in the two-dimensional planes in Figure 10. The plots show that the inferred evolution is uncertain and depends crucially on the intrinsic scatter, especially when applying a luminosity evolution for the host (see Appendix A). Assuming the lognormal prior σ_int, to mimic the assumption of the method discussed in Section 2.3, one sees that the best estimate of γ is positive, but the 95% confidence intervals extend to zero. Thus, one cannot conclude that evolution is significantly detected in our data. When relaxing the prior on σ_int, we find that the scatter is consistent with being as low as in the local samples, with a one-sided interval including zero and 1σ upper limit at around 0.5 dex. We also study the selection effect by only considering the 32 new AGNs, resulting in a higher evolutionary trend with a higher value of γ. We show the results of γ and σ_int in Table 5.

Zoom In Zoom Out Reset image size

A simple check using our prior estimate of the bias yields a similar result. The mean offset of the sample of 32 objects from the local relationship is 0.43 ± 0.06 dex, which would correspond to γ = 1.08 ± 0.15. However, after correcting for the selection bias (0.21 dex; see Section 2.3), the offset reduces to 0.22 ± 0.06 and thus γ = 0.55 ± 0.15, marginally positive but not conclusively inconsistent with the local value given the error bars and uncertainty in the correction. We also note that the distribution of offsets shown in Figure 9 displays a positive asymmetric tail. We expect the negative tail to be suppressed by our selection function, and thus it should be accounted for in our treatment. If one computes the median offset instead of the mean, the offset is almost completely consistent with the expected bias. Larger samples with different selection functions are needed to establish whether the positive tail is real or due to small sample statistics.

The local sample of inactive galaxies used in this analysis is mainly comprised of bulge-dominated galaxies, and the entire local M_* we adopted is their bulge masses. That is, in previous sections, we are comparing the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,total relations in the distant universe to the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,bulge relations locally. Considering that a significant stellar component of the high-z AGNs have a disk component (i.e., AGN hosts with fitted Sérsic index close to 1), M_*,bulge must be smaller than M_*,total. Given that the structures of our AGN hosts are similar to inactive galaxies at equivalent redshifts and stellar masses (Section 5.1), we expect their B/T ratios to follow a similar distribution that can be estimated using a single Sérsic index (Bruce et al. 2014). We can then infer the bulge stellar mass of our AGN hosts based on the inferred B/T ratios.

First, we establish a relation between the single Sérsic index and the B/T ratio using inactive galaxies from CANDELS at similar redshifts and stellar masses (Figure 11; $\mathrm{log}({M}_{* }/{M}_{\odot })\in \,[9.0,11.5]$) using the B/T measurements of Dimauro et al. (2018). In the figure, the red line indicates the average B/T ratio at a given Sérsic index. We then implement a Monte Carlo approach by randomly sampling a Gaussian Sérsic distribution for each AGN in our sample based on our measurements with 1σ errors (Table 4). Next, we randomly sample the associated B/T ratio for a given Sérsic index. To avoid the case of an unphysical faint bulge flux, we set the lower limit for the B/T ratio as 0.1. In each realization, we compare the M_*,bulge to the ${{ \mathcal M }}_{\mathrm{BH}}$ and estimate the offset and γ. We use 10,000 realizations to make sure the distribution of random samples is stable. Note that the Monte Carlo approach would take the scatter into account, including the Sérsic index and its relation with the B/T ratio. We list the resulting bulge masses M_*,bulge in Table 4.

Zoom In Zoom Out Reset image size

As expected, we find a stronger offset of ${\rm{\Delta }}\mathrm{log}{{ \mathcal M }}_{\mathrm{BH}}^{\mathrm{bulge}}\,=0.87\pm 0.07$ dex in ${{ \mathcal M }}_{\mathrm{BH}}$/M_*,bulge with respect to the local relation, which corresponds to an evolutionary trend with γ = 2.09 ± 0.30. Since selection effects are independent of B/T in our framework, we can apply our prior estimate of the correction for selection effects by removing 0.21 dex, still leaving a significant offset of 0.67 ± 0.07 dex. We show the histogram of ${{ \mathcal M }}_{\mathrm{BH}}$/M_* and γ inferred by the Monte Carlo approach in Figure 12 and illustrate the offset between bulge and black hole as a function of redshift in Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even with the uncertainty associated with our estimate of the bulge-to-total ratio, it is clear that there is a significant evolution in the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,bulge relation, considering the marginal evidence for evolution in the ${{ \mathcal M }}_{\mathrm{BH}}$–M_*,total relation and that our galaxies are not pure bulges. This finding confirms with higher fidelity the conclusion by Bennert et al. (2011b) based on optical data of a smaller sample of 11 objects in the same redshift range (see Jahnke et al. 2009; Cisternas et al. 2011; Schramm & Silverman 2013, for similar results). We list the inferred γ using the different sample combinations in this section in Table 6.

In this work, we use state-of-the-art techniques to measure the host galaxy flux, separating it from that of the point source. The fidelity of the host galaxy magnitude is high, as demonstrated by the robustness with respect to the choice of PSF, which is the primary source of uncertainty. However, we introduce some uncertainty by adopting two common simple stellar population templates to derive the rest-frame R-band luminosity and stellar mass for the full population. In principle, if we had more information on the color, we could improve on this source of uncertainty by adopting a specific stellar population model for each target. However, considering that the host magnitude in the HST/ACS band is faint (host-to-total flux ratio <30%), the individual color that we measure is insufficient to further discriminate between templates; thus, we adopt the common templates for simplicity and to facilitate reproduction of our work.

In any case, we stress that the uncertainty in the scaling relations is dominated by single-epoch black hole mass estimates (i.e., 0.4 dex), and not by uncertainties in the photometry or stellar mass estimates (i.e., typically ∼0.15 dex). Thus, the choice of template is a subdominant source of error. Note that the large uncertainties of ${{ \mathcal M }}_{\mathrm{BH}}$ are not immediately apparent from the figures, since only the AGNs that have measured ${{ \mathcal M }}_{\mathrm{BH}}$ that falls into our selection window (i.e., log(${{ \mathcal M }}_{\mathrm{BH}}$ $/{M}_{\odot })\in [7.5,8.56]$) have been selected. A posteriori, the fact that the intrinsic scatter we observe is not larger than the one in the local universe is consistent with this statement.

To compare our high-z sample to the local relation, we adopt the local measurements by Bennert et al. (2010, 2011a) and Häring & Rix (2004). In the literature (Kormendy & Ho 2013; hereafter K13; Bentz & Manne-Nicholas 2018), other analyses of the local sample are available, which could also be considered as the local anchor. However, from the point of view of our differential evolutionary measurement, we stress once again the importance of selecting local anchors that have been measured and calibrated in a similar way to the distant high-z samples. It is crucial that the black hole masses be obtained in the same of self-consistent manner, and it is also critical that the host/AGN decomposition be measured in the most similar way possible.

Therefore, for these reasons, the B10 and HR04 samples are the most appropriate for our goal, in spite of their limitations. For instance, B10 did not consider morphological features like strong bars in the decomposition. This choice might affect the inference of the bulge/disk properties in absolute terms, but it is the most similar decomposition to what is possible at high-z, where these features cannot be resolved. Also, the HR04 sample uses Jean modeling to estimate the M_* rather than using color information. This gives rise to systematic uncertainties related to the choice of the IMF, for example, which need to be kept in mind. We refer the reader to K13 for more discussion of systematic issues in local samples.

The issue of a consistent calibration of black hole masses is perhaps the most subtle and thorny. Not all of the local correlations are mutually consistent, especially when considering different galaxy properties, like luminosity, stellar mass, and velocity dispersion, and therefore a single self-consistent calibration is not possible. To illustrate the problem, we compare our local ${{ \mathcal M }}_{\mathrm{BH}}$–M_* sample to that published by K13. We find that the ${{ \mathcal M }}_{\mathrm{BH}}$–M_* measurements by K13 have a global ∼+0.3 dex offset relative to our local sample. Note that a +0.3 dex difference would almost erase the ${\rm{\Delta }}\mathrm{log}$ ${{ \mathcal M }}_{\mathrm{BH}}$ for the high-redshift sample, as described in the last section, if it were independent of the black hole mass calibration. However, the ${{ \mathcal M }}_{\mathrm{BH}}$–σ_* relations by K13 are also offset by ∼+0.3 dex with respect to those given by Woo et al. (2010), which is the baseline for our adopted black hole mass estimators. That is, if we were to adopt the K13 normalization instead, the ${{ \mathcal M }}_{\mathrm{BH}}$ of all of the AGN sample should be increased by ∼0.3 dex at the same time. As a result, adopting a different local anchor would only globally affect the level of ${{ \mathcal M }}_{\mathrm{BH}}$ to the local and high-z samples at the same time, with the result of leaving ${\rm{\Delta }}\mathrm{log}$ ${{ \mathcal M }}_{\mathrm{BH}}$ unchanged. To illustrate this point further, we adopt the K13 local ellipticals and classical bulges as the z ∼ 0 anchor and compare to the recalibrated ${{ \mathcal M }}_{\mathrm{BH}}$ of our high-z AGN sample. The inference of γ changes to 0.74 ± 0.31, which is consistent at a 1σ level with the value (i.e., γ = 1.03 ± 0.25) based on our chosen anchor.

Bentz & Manne-Nicholas (2018), given the superior resolution of their data in the local universe, decomposed their host using more components than bulge+disk (e.g., bar(s), barlens, ring(s)). This approach is impossible at high-z given current data; therefore, their sample should not be used as a local anchor, for consistency. Just for illustration, when adopting Bentz & Manne-Nicholas (2018) as the local anchor, we obtain γ = 2.12 ± 0.32 using the 32 AGNs, still indicating a positive evolution, although, as expected, the amount of evolution is different for the reasons given above.

In conclusion, whereas the uncertainty in the local anchor does not affect our measurement of evolution, it is important to keep this issue in mind in evolutionary studies and verify that the measurements are as self-consistent as possible in the local and distant universe. In order to further reduce the uncertainty, one should perform a self-consistent measurement with identical techniques and data quality in both the local and distant universe. This effort is currently in progress and, when completed, should enable us to completely eliminate this source of uncertainty (Bennert et al. 2011b, 2015; Harris et al. 2012).