The Sloan Digital Sky Survey Reverberation Mapping Project: The Black Hole Mass–Stellar Mass Relations at 0.2 ≲ z ≲ 0.8

We measure the correlation between black hole mass M BH and host stellar mass M * for a sample of 38 broad-line quasars at 0.2 ≲ z ≲ 0.8 (median redshift z med = 0.5). The black hole masses are derived from a dedicated reverberation mapping program for distant quasars, and the stellar masses are derived from two-band optical+IR Hubble Space Telescope imaging. Most of these quasars are well centered within ≲1 kpc from the host galaxy centroid, with only a few cases in merging/disturbed systems showing larger spatial offsets. Our sample spans two orders of magnitude in stellar mass (∼109–1011 M ⊙) and black hole mass (∼107–109 M ⊙) and reveals a significant correlation between the two quantities. We find a best-fit intrinsic (i.e., selection effects corrected) M BH–M *,host relation of log(MBH/M⊙)=7.01−0.33+0.23+1.74−0.64+0.64log(M*,host/1010M⊙) , with an intrinsic scatter of 0.47−0.17+0.24 dex. Decomposing our quasar hosts into bulges and disks, there is a similar M BH–M *,bulge relation with slightly larger scatter, likely caused by systematic uncertainties in the bulge–disk decomposition. The M BH–M *,host relation at z med = 0.5 is similar to that in local quiescent galaxies, with negligible evolution over the redshift range probed by our sample. With direct black hole masses from reverberation mapping and the large dynamical range of the sample, selection biases do not appear to affect our conclusions significantly. Our results, along with other samples in the literature, suggest that the locally measured black hole mass–host stellar mass relation is already in place at z ∼ 1.


INTRODUCTION
The observed scaling relations between supermassive black hole (BH) masses and the properties of their host galaxies (e.g., stellar mass and stellar velocity dispersion) in the local Universe are the foundation of modern BH−galaxy co-evolution models (Magorrian et al. 1998;Ferrarese & Merritt 2000;Gebhardt et al. 2000;Häring & Rix 2004;Gültekin et al. 2009;McConnell & Ma 2013;Kormendy & Ho 2013, and references therein).
The tight correlations suggest that active galactic nuclei (AGNs) may play important roles in regulating star formation in the host galaxies via self-regulated BH growth and feedback processes (Silk & Rees 1998;Di Matteo et al. 2005;Heckman & Best 2014).Studying BH scaling relations beyond the local Universe is a key to understanding BH and galaxy (co-)evolution over cosmic history.
The measurements of BH mass−host relations can be challenging beyond z ∼ 0.1 for several reasons.First, direct BH mass measurements based on resolved stellar/gas dynamics are difficult to obtain beyond the local Universe where the BH sphere of influence cannot be readily resolved.Reverberation mapping (RM; Blandford & McKee 1982;Peterson 2014) is the primary method of measuring BH masses for broad-line (BL) AGN beyond the local Universe, but RM is resourceintensive and only available for a small number of objects beyond z ∼ 0.1 (e.g., Bentz et al. 2013).A secondary BH mass recipe, the single-epoch (SE) virial estimator, is based on the broad-line region (BLR) radius−luminosity relation (the R − L relation) and can be easily adapted for large samples of BLAGNs at higher redshift.However, Shen & Kelly (2010) demonstrated that there is a statistical bias in SE BH masses for flux-limited samples from the uncertainties in these BH masses.In addition, the applicability of SE masses to the high-redshift and high-luminosity regime is not well understood, primarily because the local RM AGNs used to derive the R − L relation is not representative of the general quasar population (e.g., Shen et al. 2015b;Fonseca Alvarez et al. 2020), and the extrapolated R − L relations for broad Mg ii and C iv used for high-redshift BLAGNs are not as well-studied as the local R − L relation based on broad Hβ (Bentz et al. 2013).
Host-galaxy properties are also difficult to measure as the unobscured AGN (where virial BH masses are feasible) usually far outshines the host galaxy.For imaging studies, high-resolution images, such as those from the Hubble Space Telescope (HST), are often necessary to robustly decompose the quasar and host light.However, rigorous image analysis reveals that host galaxies of local AGNs (z < 0.35) often consist of complex structures, including spiral arms, tidal and merger features, in addition to the main galaxy components (bars, bulges, and disks) (Kim et al. 2017).These complex structures are extremely challenging to measure even with HST at higher redshifts.
Due to difficulties in obtaining BH mass and host properties, many studies are limited to specific samples that may introduce selection biases.Earlier studies were often restricted to the bright end of BLAGN, have small sample sizes and limited dynamical ranges in BH/host properties.Lauer et al. (2007) showed that over-massive BHs are favored in flux-limited studies due to the intrinsic scatter of the scaling relations.For a "bottom-heavy" galaxy luminosity function, there are more low-mass hosts than high-mass ones.However, more massive BHs are preferentially selected in a fluxlimited sample based on AGN luminosity, resulting in an average offset in the BH mass−host relations, and a shallower slope than the true underlying relation.Schulze & Wisotzki (2011, 2014) argued that additional selection biases could arise from the lack of knowledge in the relevant underlying distribution functions (e.g., the active fraction of AGNs, bulge properties) and their evolution with redshift.These biases can account for a large portion of, if not all, the redshift evolution reported in earlier investigations (Schulze & Wisotzki 2011;Shen et al. 2015a).
In this work, we study the BH scaling relations at 0.2 z 0.8 using the Sloan Digital Sky Survey Reverberation Mapping (SDSS-RM, Shen et al. 2015b) sample.The SDSS-RM sample has two major advantages in measuring the redshift evolution of BH−host galaxy relations: (1) the parent sample is a uniformly selected flux-limited BLAGN sample, thus the selection effects can be quantified and corrected, and (2) BH masses are available from direct RM (Shen et al. 2016;Grier et al. 2017Grier et al. , 2019;;Homayouni et al. 2020), rather than from SE, masses.We have acquired high-resolution imaging for the SDSS-RM sample with HST to measure the host galaxy color and luminosity in two bands, tracing young and old stellar populations, respectively.Our sample includes 38 sources (10 included in a pilot study in Li, J. I. et al. 2021) et al. 2013), and has sufficient statistics and dynamic range in BH mass (and stellar mass) to characterize the redshift evolution of BH scaling relations over the redshift range of 0.2 z 0.8.This paper is organized as follows.We describe our data and analysis in Section 2. The main results are presented in Section 3. We discuss our results in Section 4 and conclude in Section 5. Throughout this paper we adopt a flat ΛCDM cosmology with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 .All host-galaxy measurements refer to the stellar population only.

Sample
Our sample consists of 38 SDSS-RM quasars at 0.2 z 0.8 (median redshift z med = 0.5) with RM-based BH masses; 37 of these RM masses were based on the broad Hβ line (Shen et al. 2016;Grier et al. 2017), with one source (RM767) based on the broad Mg ii line (Shen et al. 2016).Ten sources in our HST sample were studied in a pilot program (Li, J. I. et al. 2021); 28 sources are presented in this work for the first time.Among the 44 quasars with Hβ RM BH masses in Grier et al. (2017), seven sources beyond z ∼ 0.8 were excluded from the HST programs to ensure more robust host-galaxy measurements and to avoid unknown selection biases, as the lag-detection fraction at z 0.8 is significantly lower than that at lower redshifts (e.g., see Figure 1).Figure 1 presents the redshift and luminosity distribution of our sample, and Table 1 summarizes the physical properties of these objects.

Black-Hole Masses
Reverberation mapping determines BH masses by measuring the time delay in variability between the continuum and broad emission lines.The time delay corresponds to the light travel time between the continuumemitting accretion disk and the BLR.Assuming the BLR is virialized, a BH mass can be calculated using the average time lag (τ ) and the width of the broad emission line (∆V ) via the equation: where G is the gravitational constant and f is a dimensionless factor of unity order that accounts for BLR geometry, kinematics, and inclination.The line width ∆V can be estimated from either the full-width halfmaximum (FWHM) or the line dispersion (σ line ) of the broad line measured from the mean or RMS spectra (e.g., Wang et al. 2019).
For the majority of our sources, we adopt the RM black hole masses from Grier et al. (2017) computed using a consant virial coefficient of f =4.47 based on σ line measured from the RMS spectra (equivalent to f = 1.12 when using the FWHM for ∆V ).For two of our sources, RM316 and RM519, the original σ line measurements in Grier et al. (2017) based on the first-year SDSS-RM spectroscopy are significantly overestimated; we adopt updated σ line measurements based on the 4-year SDSS-RM spectroscopy for these two objects.For RM767, Shen et al. (2016) identified a lag between the continuum and broad Mg ii line during the first-year monitoring.However, the lag significance is reduced in the more recent analysis in Homayouni et al. (2020) using 4-year light curves, as the broad Mg ii line does not display strong response to the continuum in the following years.We adopt the Shen et al. (2016) Mg ii lag for RM767, and use its σ line measured from the RMS spectrum to derive a BH mass.The BH mass uncertainties are calculated by propagating the statistical uncertainties of the lag and line width measurements, and then adding a systematic uncertainty of 0.16 dex, which is the scatter estimated from repeated RM measurements in local RM campaigns (Fausnaugh et al. 2017).However, the adopted BH mass uncertainty is still an underestimation, as it does not account for the intrinsic scatter in the virial coefficient for individual systems, which could lead to additional BH mass uncertainties of as large as ∼ 0.3 dex.The BH masses are tabulated in Table 1 (with updates from earlier work indicated by an asterisk).

HST Imaging Analysis
The HST observations for the 28 new objects were conducted between 2019 December 23 and 2021 June 09  Davis, Davis et al. 2007).We processed the individual UVIS exposures from this earlier program following the same procedures for our HST program, and the final imaging sampling is 0 .05pixel −1 (∼ 0.3 kpc at z = 0.5).We use a field star in the same field of view as the PSF model for the UVIS images of RM177.
We then follow the procedures in Li, J. I. et al. ( 2021) and perform 2D image decomposition to separate the quasar and host light using GALFIT (Peng et al. 2010), but with two modifications.First, we allow each system to be fitted by a PSF+disk model (Sérsic index n = 1), in addition to the PSF+bulge model (n = 4) and PSF+bulge+disk model.Upon our analysis with the full sample, we identified several sources whose hosts best fitted by an exponential disk, rather than a bulge or bulge+disk model as used in Li, J. I. et al. (2021).Second, we revise our model-selection criteria using reduced-χ 2 calculated from a small region surrounding the target.By default, GALFIT's reduced-χ 2 is calculated from the entire image analysis area (e.g., roughly 10 ×10 ) where > 60% of all the pixels are background, so the reduced-χ 2 can change based on the chosen image size, and is largely determined by the accuracy of the background estimation.To assess the quality of the fit, we calculate rχ 2 e , the reduced-χ 2 within the best-fit ellipse at 3σ sky background (estimated by GALFIT) around the source.If ∆(rχ 2 e ) > 5 (threshold chosen by visual inspection of the data) between the PSF+bulge+disk model and the 2-component models (PSF+bulge or PSF+disk), we consider there is strong evidence that the additional component is necessary and adopt the three-component model; otherwise, we select the two-component model (PSF+bulge or PSF+disk) with the smaller rχ 2 e as the best-fit model.To briefly summarize our fitting procedure, we first fit the IR images with three different models: PSF+bulge (Sérsic index n = 4), PSF+disk (n = 1), and PSF+bulge+disk.The fit is considered successful when the best-fit parameters are within reasonable ranges (i.e., the effective radius of the Sérsic component R e > 1 pixel, axis ratio q > 0.01), which is to prevent introducing additional components fitting for mismatched PSF or other small-scale features.While the galaxy may not be a perfect bulge or disk, Kim et al. (2008) showed that fixing the Sérsic index results in more accurate flux recovery during host decomposition when the host galaxy is faint.We use ∆(rχ 2 e ) to select the bestfit model from the successful PSF+bulge, PSF+disk, and PSF+bulge+disk models.
In addition to the quasar+host, we fit additional PSF and/or Sérsic models for nearby objects to ensure the host decomposition and sky background estimation are not strongly affected by nearby objects (e.g., see RM033, RM101, RM694, RM776, etc, for examples).
We visually inspect all the GALFIT images and manually adjust the GALFIT models only when necessary.Upon visual inspection, the background in RM776 is high due to a nearby bright object, and adding another component improves the fitting of its surfacebrightness profile significantly, so we adopt a threecomponent model for RM776.RM775 and RM790 display extended truncated ring features in the residual images of the PSF+bulge+disk model, so a fourth component (an inner-truncated disk) was added to ensure robust flux recovery for the host.The truncated disk in RM775 is also fitted with Fourier modes to account for the irregular ellipsoid shape.However, we only include the main disk component in the PSF+bulge+disk model, and not the truncated disk, for estimating the final photometry for the disks.Finally, we fit the flux of each component in the UVIS images by fixing the shape and structural parameters (Sérsic index, effective radius, ellipticity, and position angle) to the best-fit model in the IR images.For the sources that preferred the threecomponent model in the IR image, we check if the three components in the UVIS image converge on similar relative positions as in the IR image, which ensures the model is fitting the same physical structures in the two bands.The bulge and disk components of two sources, RM267 and RM316, failed to converge at similar central Li et al. positions, so the two-component model (PSF+bulge) is adopted instead.Fig. 2  During our analysis of the full HST sample, we discovered an error in our GALFIT analysis in the pilot study (Li, J. I. et al. 2021).The ncombine parameter was input incorrectly, which caused the sigma image produced by GALFIT to be overestimated by a factor of ∼ 4 in areas dominated by emission (see GALFIT user manual, Equation 33).The error mainly affects the estimation of χ 2 , but does not change the fitting results, i.e., all fitted parameters are consistent with the results with the correct sigma images within the uncertainties.We include updated measurements for the 10 objects in the pilot study in Table 2.
GALFIT only accounts for statistical uncertainties between the data and the model, and does not take into account PSF mismatches or complex spatial structures.There are three major sources of flux uncertainties: (1) the temporal variability of the HST PSF (derived from the difference between the dedicated PSF observation and field stars in science observations, ∼ 0.07 mag in UVIS and ∼ 0.03 mag in IR), (2) the deviation between the GALFIT model and the image (∼ 0.02 mag in UVIS and ∼ 0.005 mag in IR), and (3) fixing the Sérsic index (∼ 0.05 mag for PSF and ∼ 0.2 mag for the host/bulge/disk).We combine these flux uncertainties and adopt typical values of 0.1 and 0.25 mag as the final uncertainties for the PSF and galaxy (bulge, disk, or galaxy) flux measurements in all bands, respectively.These final uncertainties are consistent with those in our pilot study and similar observations and simulations in the literature (e.g., Kim et al. 2008;Jahnke et al. 2009;Park et al. 2015;Bentz & Manne-Nicholas 2018).See Li, J. I. et al. (2021) for additional technical details on the flux uncertainty budget.

Host Galaxy/Bulge Masses
Following the approach in Li, J. I. et al. (2021), we convert the UVIS/IR photometry to rest-frame B and I/R band and estimate the host/bulge stellar masses with the color-M * /L relations (CMLR) from Into & Portinari (2013) and CIGALE (Boquien et al. 2019).First, we correct for Galactic extinction using the recalibrated Schlegel et al. (1998) dust map and reddening from Schlafly & Finkbeiner (2011).We then fit the extinction-corrected HST photometry with CIGALE to derive k-corrections and color transformations between the HST filters and the Johnson−Cousins filters.CIGALE is a spectral energy distribution (SED) fitting code that can model galaxy and AGN emission from multiwavelength photometry.We set up a simple CIGALE model that includes basic stellar population synthesis models (Maraston 2005), a initial mass function (Kroupa 2001), a dust attenuation model (Calzetti et al. 2000;Leitherer et al. 2002), and a delayed star formation history with optional starburst.We do not include the AGN model for modeling the quasar-subtracted photometry.The UVIS filters are converted to B-band magnitudes, and the IR filters are converted to I-band and R-band filters depending on the source redshift (F110W to I-band at z < 0.4 and R-band z > 0.4; F140W to I-band at z < 0.7 and R-band at z > 0.7).
We estimate the host and bulge stellar masses with the CMLR for dusty galaxy models from Into & Portinari (2013) using the rest-frame photometry and their uncertainties.CIGALE fits provide the k-corrected photometry, from which we estimate stellar mass with the CMLR relation, and a stellar mass from the best-fit SED model.For RM177 (with two UVIS bands from a separate HST program), we include both the F606W and F814W bands for the CIGALE fitting but only use the F606W band (rest-frame R-band) for the CMLR stellar mass estimation.The CMLR stellar-mass uncertainties are propagated directly from the photometry uncertainties, and the CIGALE stellar-mass uncertainties are estimated from the SED modeling.Both CMLR and CIGALE uncertainties are consistently around 0.3 dex, which is typical for stellar mass estimation from twoband photometry.
The final Galactic-extinction corrected, k-corrected, band-converted magnitudes, and the host/bulge stellar masses are tabulated in Table 3.We adopt the CMLR stellar masses as our nominal host/bulge stellar masses.The best-fit stellar masses from CIGALE are also reported for comparison, which are generally consistent with those estimated from the CMLR.The only exception is RM265.The color derived for RM265 from CIGALE is unusually red, which led to a large, likely unphysical, host stellar mass (> 10 12 M ) using the CMLR.However, the typical B − R color is roughly 0.3 < (B − R) < 1.3, derived from all galaxy types in the Kinney-Calzetti Spectral Atlas (Calzetti et al. 1994;Kinney et al. 1996).If we assume a red color of 1.3 and adopt the R-band luminosity for the CMLR, the host stellar mass for RM265 is log(M * ) = 11.23, which is consistent with the stellar mass derived from CIGALE.We show both the CMLR and CIGALE masses for RM265 in Figures 5, 6, and 10, and use the more physical CIGALE mass (for RM265 only) when fitting the BH scaling relations and their redshift evolution in our analysis.

Host Properties
At z > 0.2, it becomes challenging to perform bulge/disk decomposition due to limited spatial resolution, even with HST.Our GALFIT analysis shows 16 (out of 38) quasars are best-fitted by the PSF+bulge model, nine quasars are best-fitted by the PSF+disk model, and 13 quasars are decomposed into PSF+bulge+disk models.In addition, 26 hosts are bulge-dominated, i.e., M * ,bulge > M * ,disk , and 12 hosts are disk-dominated.A best-fit profile of n = 4 (n = 1) in our analysis does not necessarily mean the host galaxy is an elliptical (spiral) galaxy; the Sérsic index is fixed to n = 1 or n = 4 to ensure the quasar/host decomposition is robust and not to provide rigorous classifications of host morphology.In fact, the majority of local elliptical galaxies are not well-described by single Sérsic components (e.g., Huang et al. 2013), and exponential profiles do not always indicate the presence of disks.
The structural parameters (ellipticity, Sérsic index, effective radius) of the bulge/disk-dominated sources in our sample are broadly consistent with the statistical distributions from ∼ 2500, i-mag<22 SDSS quasar hosts observed by the Hyper Suprime-Cam (HSC) on the Subaru telescope (Li et al. 2021b).When we allow the Sérsic index to vary in the GALFIT fitting, the median (minimum, maximum) Sérsic index of our two-component model is 2.0 (0.6/7.0), similar to the distribution in Li et al. (2021b).There are more disk-like (n < 2) hosts in the SDSS-HSC sample, but roughly equal numbers of bulge-like and disk-like hosts (Sérsic indices above and below 2) in our sample.The size and ellipticity of our quasar hosts are also similar to the SDSS-HSC sample.The median (16%, 84% percentiles) effective radius is 0 .68(0 .35/0 .92),and the median (16%/84% per- Note-r is the effective radius of the Sérsic component, n is the Sérsic index, q is the ratio between the semi-minor axis and the semi-major axis, and P.A. is the position angle at the semi-major axis in degrees.The reduced χ 2 is calculated from the image residual, as reported by GALFIT.Magnitudes are reported in ST magnitude (mag which is the default output from GALFIT.For RM775 and RM790, we include the best-fit parameters for the truncated disks in the UVIS and IR images: the magnitudes are the surface brightness at the break radius (mag/arcsec 2 ), and the best-fit parameters for the truncated radial profiles are listed in the order of the 1% flux radius (softening length, in arcseconds), 99% flux radius (break radius, in arcseconds), q, and P.A..The truncated disk in RM775 is fitted with Fourier modes in both the disks and truncated radial profiles, and the best-fit Fourier amplitudes (first row) and phase angles (second row) are listed in the order of Fourier mode 1, 3, 4, 5, 6.No extinction corrections are made for these magnitudes.The uncertainties of the GALFIT results are discussed in Section 2.
We also examine the offset between the quasar position and the host centroid in the IR images, where the centroid of the host galaxy is better constrained than in the UVIS band.Off-centered AGN/quasars may indicate on-going galaxy mergers or recoiling SMBHs from binary SMBH coalescence (Loeb 2007;Comerford & Greene 2014).Figure 3 shows that most (34/38) of the quasars are located within < 1 kpc of the host galaxy center.The four sources with significant offsets (> 1 kpc, RM265, RM267, RM634, RM645; see the images and GALFIT models in full Figure 2 figure set online) show signs of galaxy interaction or mergers, which would complicate the centroid measurements of the host galaxy.These results suggest z < 1 quasars are well centered within ∼ 1 kpc of the host centroid, consistent with the findings using alternative approaches (Shen et al. 2019).
While studies of local AGNs have demonstrated that BH properties mainly correlate with the bulge and not the entire host (e.g., Kormendy & Ho 2013), studies at higher redshift are often limited to the BH−host relations when bulge/disk decomposition is difficult or impossible (e.g., Jahnke et al. 2009;Merloni et al. 2010).In this work, we present both the BH−bulge and the BH−host relations in our sample, where M * ,bulge and M * ,host refer to the bulge-only and total host stellar mass, respectively.We include all sources in the M BH − M * ,host relation and exclude the disk-only (PSF+disk) objects in the M BH −M * ,bulge relation.When comparing with earlier work, we examine their bulge/disk decomposition assumptions and place the comparison on an equal footing, i.e., including bulge-dominated or bulge/disk decomposed sources only in the M BH −M * ,bulge relation, and including all sources in the M BH − M * ,host relation.

Comparison with earlier work
Figure 4 shows the M BH −M * ,host and M BH −M * ,bulge relations of our sample and several local and higherredshift samples.High-resolution HST imaging has been used to investigate the AGN host galaxies at z > 0.2 (e.g., Jahnke et al. 2009;Bennert et al. 2011;Bentz   (2020) that estimated M * ,bulge by assigning bulge/total ratios depending on the Sérsic indices of the host profile.At even higher redshift (e.g., z 1.5), host stellar masses can be obtained by SED fitting (e.g., Merloni et al. 2010;Dong & Wu 2016;Suh et al. 2020) or imaging analysis of lensed quasars (Peng et al. 2006b;Ding et al. 2021).SED fitting with wide wavelength coverage can provide better color information for estimating the stellar masses (compared to using only two HST bands), and large samples can be studied simultaneously in multiwavelength fields.However, it is impossible to distinguish between the bulge and disk components through SED fitting.M * ,host can also be measured from the reconstructed images of strongly lensed quasars up to z ∼ 3. Our sample is generally consistent with the M BH − M * ,bulge and M BH − M * ,host relations of these intermediate-to-high redshift samples.
Our sample is the only uniformly-selected (i.e., selected based on a flux limit) AGN sample with RMbased BH masses to study the BH scaling relations beyond the local Universe (z > 0.1).The RM masses in Grier et al. (2017)   mendy & Ho (2013) using the virial factor from Woo et al. (2015).BH masses in all the comparison samples, except for Bentz & Manne-Nicholas (2018), are derived from the SE method, which is less reliable than RM masses.The SE method relies on a "tight" R − L relation to estimate BLR sizes based on quasar luminosities; however, recent studies (Du et al. 2016;Fonseca Alvarez et al. 2020) have shown the local R − L relation is biased towards the local AGN sample and could be overestimating SE BH masses by as much as ∼ 0.3 dex when applying to the general quasar population.Finally, the SE method is calibrated to local quiescent galaxies or local RM AGNs, and different virial factors may be used for different samples or broad-line species.When comparing to samples from the literature, we rescale all M BH values using f = 4.47 as in Grier et al. (2017), even for the SE masses.We also compare our results with local baseline samples from the literature, including the quiescent galaxies (mainly ellipticals, Kormendy & Ho 2013), active galaxies (Bennert et al. 2021), and RM AGNs (Bentz & Manne-Nicholas 2018).The M BH − M * ,bulge and M BH − M * ,host relations of the three local samples and our best-fit relations are consistent in slope and intercepts within uncertainties.However, Reines & Volonteri (2015) find AGN hosts follow a similar slope as lo-cal quiescent galaxies but are an order of magnitude lower in normalization for the M BH − M * ,host relation.They suggested the difference in normalization may be due to AGN activity or galaxy morphology (which is also shown in Greene et al. 2020).These results may appear contradictory at first glance; however, it is difficult to provide a straightforward comparison since these studies adopt different stellar-mass estimation methods.Bentz & Manne-Nicholas (2018) 2021) performed a systematic analysis on the evolution of M BH − M * relations in cosmological simulations of Illustris, TNG 100, TNG 300, Horizon-AGN, EAGLE and SIMBA.They find that the median/mean M BH − M * relations at 0 < z < 1 are in general agreement with observational data and there is little evolution with redshift.The observed tight correlation between BH accretion rate and star formation at 0.5 < z < 3 indicates the growth of BH and host galaxies are in sync, and the BH scaling relations should not have strong redshift dependence (Yang et al. 2019).However, the scatter in M BH − M * relations differs in these simulations, which mainly depends on the implemented subgrid physics in the simulations, e.g., the strength and efficiency of supernova and AGN feedback.Our sample consists of 38 sources spanning more than two orders of magnitude in M BH and M * ,host/bulge , which is sufficient for statistical analysis.The Pearson correlation coefficient r between M BH and M * ,host/bulge are r ∼ 0.5 with low p-values (<0.05), suggesting the BH and galaxy/bulge masses are positively correlated and the correlation is statistically significant.We use the LINMIX ERR algorithm (Kelly 2007) to perform linear regression fitting on the M BH −M * ,bulge and M BH −M * ,host relations.LINMIX ERR is a Bayesian fitting algorithm that accounts for uncertainties in both axes and intrinsic scatter in the relations.We fit for the equation: and tabulate the best-fit parameters in Table 4.The regression fits to the observed sample do not account for selection effects.In §4.2 we use a more robust fitting code to constrain the intrinsic BH-host scaling relations, and we include the bias-corrected best-fit parameters in Table 4..However, given the large dynamic   −0.17 range of our sample, selection effects do not appear to impact the results significantly, as we will show in §4.2.

Comparison with spectral decomposition
We compare the host-light fraction and total stellar mass derived from HST imaging decomposition in this work and the spectral decompositions in Shen et al. (2015a) and Matsuoka et al. (2015).Both of these earlier studies measured host-galaxy properties using the high-S/N coadded spectra from the first-year SDSS-RM spectra (Shen et al. 2015b).Shen et al. (2015a) used a principal component analysis method to decompose coadded spectra into quasar and galaxy spectra, and measured host-galaxy properties directly from the galaxy spectra, including stellar velocity dispersion and host-free AGN luminosity.Matsuoka et al. (2015) performed spectral decomposition on the coadded spectra using models of AGN and galaxy spectra, and measured host galaxy properties by fitting the decomposed galaxy spectra with stellar population models.
To compare the host-light fraction (f * , the fractional contribution of the host stellar component to the total flux), we calculate the host fraction in the HST imaging using the decomposed GALFIT models within the 2 diameter spectral aperture, and the host fraction in spectral decomposition by computing the expected flux in the F606W and F814W bandpass in the decomposed spectra from Shen et al. (2015a).Figure 6 (left panel) reveals that the host fraction from image decomposition is systematically higher than that derived from spectral decomposition, similar to our finding in the pilot study (Li, J. I. et al. 2021) and in Yue et al. (2018).Figure 6 (right panel) compares the host stellar mass derived from this work and from Matsuoka et al. (2015).Our stellar mass is systematically smaller by ∼0.5 dex.The cause of the stellar mass offset is currently unclear, but it might be partially due to different choices of initial mass functions and stellar population models in Matsuoka et al. (2015) and CIGALE, or the fiber-loss correction applied in Matsuoka et al. (2015), which assumes the mass-to-luminosity ratio in the central region (within the 2 -diameter aperture) represents that for the entire galaxy.ple due to intrinsic scatter in BH-host scaling relations (Lauer et al. 2007), we perform a forward-modeling simulation following the procedures in Shen et al. (2015a).
We first simulate a parent quasar sample following the local M * distribution from Bernardi et al. (2010) and the M BH − M * ,bulge relation from Kormendy & Ho (2013), with an intrinsic scatter of 0.29 dex.Using the true M BH , we assign a quasar bolometric luminosity by assuming a lognormal Eddington ratio distribution (λ ≡ L bol /L Edd ) and the Eddington luminosity is L Edd = 1.26 × 10 38 (M BH /M ) erg s −1 .We choose a mean Eddington ratio of logλ = −1 and a scatter of 0.3 dex (Shen et al. 2008;Shen & Kelly 2012).We include measurement uncertainties of 0.35 dex for M * and 0.2 dex for M BH to mimic the uncertainty levels of our CMLR-based M * and RM M BH measurements.Finally, for 100 bootstrap iterations, we randomly draw 38 sources to perform LINMIX ERR fitting with at different i-mag < 16, 19, and 22 (similar to the flux limit of the SDSS-RM sample).
Figure 7 shows how the flux limit biases the observed scaling relations.When the flux limit increases, overmassive BHs are preferentially selected, the slope of the best-fit relation becomes shallower, and the normalization increases.The best-fit intrinsic scatter remains roughly the same in our simulations.However, our simulation does not include the outlier population with under-massive BHs seen in observations (e.g., see Figure 4).Missing the outlier population could lead to an underestimation of the intrinsic scatter and selection bias based on the flux limit, since under-massive BHs are less likely to be selected in flux-limited surveys.Given the relatively faint flux limit of our SDSS-RM sample, selection biases do not play an important role in the measured M BH − M * relations (see the right panel of Figure 7), and we found similar relations at z med = 0.5 as the local relations.Our results are consistent with other studies that properly account for selection biases (e.g., Sexton et al. 2019;Suh et al. 2020;Li et al. 2021a).We note that the selection of our sample also depends on successful RM lag measurements, which may depend on BH mass and Eddington ratio, etc.However, since the lag-detection fraction in the Grier et al. (2017) sample is nearly uniform up to z ∼ 0.8, we assume the sample selection is not strongly affected by additional selection biases based on the quasar properties.

Quantifying the Selection Effects
To quantitatively account for the underlying galaxy properties (i.e., the galaxy mass function) and selection effects, we follow the framework of Kelly (2007) to perform a Markov chain Monte Carlo (MCMC) fitting for the intrinsic scaling relations and scatters.Specifically, we wrote a MCMC fitting code based on the Metropolis-Hastings algorithm (Metropolis et al. 1953  1970) and the statistical derivation in Section 4.1 and 5.1 of Kelly (2007), which allows us to account for the selection effect based on the dependent variable.Here, we briefly summarize the fitting algorithm and parameter setup and refer the readers to Kelly (2007) for the full, detailed mathematical derivation.Our custom fitting code is available via ftp://quasar.astro.illinois.edu/public/sdssrm/paper data/Li 2023 HST host.
As shown in Kelly (2007), when the intrinsic scatter and the uncertainties are comparable to the dynamical range of the data, the best-fit slope becomes shallower when the underlying distribution is not considered as a prior in the fitting procedure.One solution is to incorporate empirical distributions into the likelihood function, e.g., the observed local stellar mass function Φ(x) from literature (e.g., Bernardi et al. 2010).However, due to the limited dynamical range in M * and the small sample size of our data, the local galaxy mass function is not a good prior for our sample.Alternatively, Kelly (2007) suggests using a series of Gaussian functions to model the underlying distribution, which provides a flexible and empirical solution even when the underlying distribution is unknown.This is the method implemented in the original LINMIX ERR fitting algorithm, which we continue to adopt in our MCMC fitting for consistency.
When the sample selection is based on the dependent variable (i.e., M BH ), the posterior distribution and likelihood depend on an additional term (P (I = 1|θ), where θ are the model parameters) that describes the likelihood of including each data point in the observed sample based on the model parameters (for more details, see Section 5.1 in Kelly 2007).Following the same procedure as described in Paragraph 1 of this section, we estimate the expected i-band magnitude by assigning a random Eddington ratio and an Eddington luminosity based on the redshift and a range of "true" M BH for each data point.The probability of including a data point is 1 if i-mag< 21.7 and 0 if i-mag> 21.7.Finally, we calculate P (I = 1|θ) by integrating the probability of including each data point and their likelihood over a range of "true" M BH and M * given the model parameters.We adopt uninformative, flat priors for all parameters (5 < a < 10, 0 < b < 3, and 0.0001 < σ 2 < 1) and minimize the product of the likelihood and prior to compute the posterior distribution of the parameters a, b, and σ.
Figure 8 and 9 present the posterior distribution of a, b, and σ, and the best-fit values are tabulated in Table 4.As expected, the slope of the intrinsic scaling relations becomes steeper, and the normalization decreases, after correcting for the selection biases.The best-fit pa-

Evolution of BH Scaling Relations
Earlier works on M BH − M * ,bulge and M BH − M * ,host relations found that the average BH-to-host galaxy mass ratio evolves positively with redshift (Peng et al. 2006a,b;Merloni et al. 2010;Bernardi et al. 2010).However, selection biases and measurement uncertainties could yield false positives of the evolution.For example, Jahnke et al. (2009) reported that there is no evidence of evolution when they carefully choose their sample to avoid selection biases.Similarly, Suh et al. (2020) found the redshift evolution seen in Merloni et al. (2010) can be explained by the Lauer bias (Lauer et al. 2007), and there is no trend of evolution in their Xray-selected, lower luminosity sample.Using over 500 uniformly-selected 0.2 < z < 0.8 SDSS quasars, Li et al. (2021a) found a redshift evolution of the offset in the M BH −M * ,host relation that is within ±0.2 dex from zero, consistent with no significant evolution since z ∼ 0.8.
Figure 10 presents the deviation of M BH from the Kormendy & Ho (2013) M BH − M * ,bulge relation of our sample.
We fit the deviation as a function of log(z ), and the slopes and intercepts are consistent with zero, suggesting there is no redshift evolution from the local relations.The median (16th/84th percentiles) black hole/host galaxy mass ratios are M BH /M * ,bulge = 0.0037(0.0007/0.0088)and M BH /M * ,host = 0.0030(0.0007/0.0076),within 1σ un- certainty of the local value of M BH /M * ,bulge ∼ 0.005 (Kormendy & Ho 2013).
Our results agree with recent observational studies that there is limited evolution in the M BH − M * relation (Suh et al. 2020;Li et al. 2021a).The limited redshift evolution in the M BH −M * relations can be explained by the tight correlation between the BH accretion rates and star formation rates found in bulge-dominated galaxies at z = 0.5 − 3 (Yang et al. 2019), which suggests the growth of SMBH and their host galaxies and bulges are in sync since z ∼ 3.

BH−Bulge versus BH−Host Relations
In the local Universe, SMBH masses are tightly correlated with the properties of classical bulges, but not with disks or total mass of the host galaxy.However, the intrinsic scatter in the M BH −M * ,host relation of our sample is slightly smaller than that in the M BH − M * ,bulge relation.There are a few extra sources of uncertainties for the bulge-mass estimate than for the total mass estimate, which will contribute to the intrinsic scatter.First of all, bulge-disk decomposition could add significant uncertainties to the bulge mass.Most of our host galaxies are far less luminous than the quasar, and compact hosts could be near the limit of HST imaging resolution.Bulge-disk decomposition is reliable when both the bulge and disk are sufficiently bright (compared to the central quasar) and there is a distinct difference in their effective radii.Furthermore, we cannot distinguish classical bulges from pseudo-bulges or other bulge-like structures from surface-brightness decomposition, nor could we model complex bar, spiral, Li et al. and merger structures in the bulge-disk decomposition, which increases the uncertainties in bulge identification and mass estimation.Gao & Ho (2017) found that rigorous modeling of bars and innermost structures (e.g., rings and disk breaks near the bulge) is crucial to recovering bulge properties, while the modeling of spiral arms and extended disks have negligible effects.We note that some host galaxies in our sample show clear evidence of bars (e.g., RM320, RM634, etc), which are modeled as bulges (n = 4) or disks (n = 1) in our analysis, without additional bar structures.Moreover, some host galaxies show clear spiral arm features (e.g., RM371, RM772, etc), indicating the presence of disks, but are modeled as "bulges"(n = 4).Previous works (e.g., Zhao et al. 2021;Greene et al. 2020) showed latetype quasar hosts preferentially scatter below early-type hosts in the M BH −M * ,bulge relation, which is not seen in our data, suggesting our bulge-disk decomposition is not as reliable in measuring galaxy morphology.A detailed simulation of bulge-disk decomposition for AGNs with similar host and quasar properties (e.g., AGN/host flux ratio, host effective radius, Sérsic indices, and complex structures) is needed to provide quantitative uncertainty estimation, which is beyond the scope of this work.
Another possible source of uncertainty is the CMLR estimation for bulges.Recent studies have reported that compact regions around the SMBH may have denser interstellar medium, boosted star formation, and complex stellar populations (e.g., Ni et al. 2019;Kim & Ho 2019;Zhuang & Ho 2020;Shangguan et al. 2020;Yesuf & Ho 2020;Molina et al. 2021).Two-band color and the use of empirical M/L relation may not be sufficient to produce reliable estimates for the bulge stellar mass.

CONCLUSIONS
We present the M BH − M * ,bulge and M BH − M * ,host relations of 38 sources with RM-based BH masses (Grier et al. 2017) and 0.2 z 0.8 (median redshift z med = 0.5).Our sample is the first uniformly-selected sample with RM-based BH masses at z > 0.3 for studying BHhost relations, and covers two orders of magnitude in BH mass and host stellar mass.The reliable RM-based BH masses and host mass estimates from HST imaging decomposition, combined with the large sample size and dynamic range in mass, allow one to alleviate selection biases in studying the potential evolution of the BH-host scaling relations.Our scaling relations are consistent with those for local AGNs, quiescent galaxies, and other high-redshift samples, with negligible redshift evolution up to z 1.As shown in Table 4, the best-fitting intrinsic M BH − M * ,host relation is: log(M BH /M ) = 7.01 +0.23  −0.33 + 1.74 +0.64 −0.64 log(M * ,host /10 10 M ) after correcting for the underlying sample distribution and selection effects.We estimate an intrinsic scatter of 0.59 +0.23  −0.21 dex and 0.47 +0.24  −0.17 dex in the M BH − M * ,bulge and M BH − M * ,host relations, respectively, which is again consistent with the local BH scaling relations.
With our approved Cycle 1 JWST proposal (GO-2057, PI: Shen), we will continue to explore BH−host relations and their redshift evolution up to z ∼ 2 using quasars with direct RM-based BH masses (Grier et al. 2019).

Figure 2 .
Figure 2. Examples of surface-brightness decomposition of three quasars with GALFIT, from top to bottom are sources that are best-fitted by a PSF+bulge+disk, PSF+disk, and PSF+bulge model, respectively.The left panels are the surface-brightness profiles of the data (black dots), the model (grey solid line) and each modeled component (red solid lines for PSFs, orange dotted-dash lines for bulges (n=4), blue dash lines for exponential disks (n=1), and purple dotted lines for truncated rings in RM775 and RM790).The radial profiles are directly-measured from the GALFIT decomposed models and the HST images with isophote fitting.The leftmost, bottom sub-panel for each object is the residual of the surface-brightness profile, with the rms along the isophote elliptical plotted in grey.The three images on the right are (from left to right) the HST image, the GALFIT model, and the residual.The blue ellipse in the HST image (IR only) encloses the area above 3σ sky background in the best-fitting model.The residual images display the 1 st to 99 th percentiles (with linear stretch) of the residual values to provide better visual contrast.The reduced χ 2 of the model is labeled in the lower right corner of each residual image.The full figure set is available online.

Figure 3 .
Figure 3. Projected physical offset between the quasar position and host centroid in the IR images as a function of the host effective radius.The uncertainties are estimated from GALFIT (quasar position and effective radius) and the quasarsubtracted host images (centroid position of the host).The error bars are inflated by a factor of 10 for clarity, as the GALFIT uncertainties are small and likely underestimated (median uncertainties are ∼0.02kpc for the quasar/host separation and ∼0.005 kpc for the effective radius).

Figure 5 .
Figure 5. BH mass as functions of bulge stellar mass (left) and total stellar mass (right) of our sample (black circles for bulge-dominated sources, gray squares for disk-dominated sources, and the open gray square for the CIGALE stellar mass of RM265).The blue dashed lines (and the gray shaded areas) are the best-fit relation (1σ range) of our sample.The red solid lines are the Kormendy & Ho (2013) local MBH − M * ,bulge relation.The Pearson-r coefficient, p-value, and intrinsic scatter of the relations are labeled in the top-left corner of each panel.ferentevolutionary paths or feedback mechanisms for establishing the BH scaling relations.Local studies of BH scaling relations also found that late-type galaxies follow a similar slope in the BH scaling relations, but at a lower normalization than early-type galaxies(Reines & Volonteri 2015;Greene et al. 2020;Zhao et al. 2021).In addition, there is no strong evolution in the M BH − M * ,bulge relation up to z ∼ 1 in the Illustris simulation.Mutlu-Pakdil et al. (2018) studied the M BH −M * ,host relations in the Illustris simulation to provide a better comparison for high redshift observations, and reported that the M BH −M * ,bulge and M BH −M * ,host relations are generally consistent with each other up to z ∼ 1.Volonteri et al. (2016) studied the M BH − M * ,host and M BH −M * ,bulge relations in the Horizon-AGN simulation.By identifying classical bulges in their simulation through kinematics and bulge/disk decomposition, they reproduced the tight M BH − M * ,bulge relation of classic bulges fromKormendy & Ho (2013).Other simulations, e.g., MassiveBlack-II(Khandai et al. 2015), generally produce similar trends in BH scaling relations at z < 1.Habouzit et al. (2021) performed a systematic analysis on the evolution of M BH − M * relations in cosmological simulations of Illustris, TNG 100, TNG 300, Horizon-AGN, EAGLE and SIMBA.They find that the median/mean M BH − M * relations at 0 < z < 1 are in general agreement with observational data and there is little evolution with redshift.The observed tight correlation between BH accretion rate and star formation at 0.5 < z < 3 indicates the growth of BH and host galaxies

3. 3 .
M BH − M * ,host and M BH − M * ,bulge Relations of our sample

Figure 6 .
Figure 6.Left: Comparison of the derived host light fraction from this work and Shen et al. (2015a), both measured in the UVIS (F606W or F814W) bandpass.The black dashed line shows the 1:1 ratio line.Right: Comparison of the total stellar masses derived from this work and Matsuoka et al. (2015) for overlapping objects.

Figure 7 .
Figure 7. Results of the Lauer et al. (2007) bias simulation described in Section 4.1.Left: grey contours represent the underlying quasar sample, and the grey points are one realization of the mock sample at the i = 22 flux limit, after being perturbed with measurement uncertainties.The black solid line is the unbiased Kormendy & Ho (2013) relation.The colored lines indicate the results of the mock sample at different flux limits.Right: the best-fit MBH − M * ,bulge parameters at each flux threshold and the best-fit MBH − M * ,host parameters of our sample (plotted in red at i = 21.7).The dotted lines show the parameters and intrinsic scatter from the Kormendy & Ho (2013) relation.

Figure 8 .
Figure 8. MCMC fitting for the MBH − M * ,bulge relation when considering the selection bias.Left: the MBH − M * ,bulge relation with the 1σ range drawn from the posterior (black dashed line and shaded area).The blue dashed lines and shaded area show the original best-fit relations from LINMIX ERR (previously shown in Figure 5), and the red solid line is the local MBH − M * ,bulge from Kormendy & Ho (2013).Right: the posterior distribution of the parameters a, b, and σ.The best-fit parameters from LINMIX ERR and the MCMC fitting are indicated by the blue and white squares, respectively.The MCMC fitting recovers a similar BH-galaxy relation to the original linear regression, indicating that the selection effects are minimal.

Figure 9 .
Figure 9. Same format as Figure 8 but for the MBH − M * ,host relation.rameters are within uncertainties as the LINMIX ERR fit (without considering selection bias) discussed in Section 3, demonstrating that our results are not strongly affected by selection biases.Intrinsic scatter of the M BH − M * ,bulge and M BH − M * ,host relations is an important indicator for BH−galaxy co-evolution, as it might be related to the galaxy/AGN properties and their evolutionary path.The local samples of Kormendy & Ho (2013) and Bennert et al. (2021) only include classical bulges and pseudo-bulges and have a smaller intrinsic scatter of 0.28 − 0.39 dex.However, when including all morphological types and active/inactive galaxies, the intrinsic scatter increases to ∼0.5 dex (Reines & Volonteri 2015; Bentz & Manne-Nicholas 2018) for the M BH − M * ,bulge relation, and becomes even slightly larger for the M BH − M * ,host relation.In addition, the BH accretion rate is found to be correlated with other host properties, e.g., compactness of the central ∼1 kpc region (Ni et al. 2019, 2021), which can introduce additional scatter in the BH scaling relations.For our quasar sample, the intrinsic scatter of the M BH − M * ,host and M BH − M * ,bulge relations are 0.47 +0.24 −0.17 dex and 0.59 +0.23 −0.21 dex, respectively, after accounting for the selection effects, which are comparable to the scatter in the local relations.The intrinsic scatter of our M BH − M * ,host relation is slightly smaller ( 0.5σ of difference) than our M BH − M * ,bulge relation, which we will further discuss in Section 4.4.Because we have neglected the systematic uncertainty in our RM BH masses due to the scatter in individual virial coeffi-

Figure 10 .
Figure 10.Evolution of ∆log(MBH) with redshift, with baselines adopted from the best-fit relations of MBH − M Bulge and MBH − L Bulge from the Kormendy & Ho (2013) sample.Our work is the only sample with RM-based BH masses beyond z > 0.3.Vertical error bars are from uncertainties in BH mass only.
and the China Manned Space Project (CMS-CSST-2021-A04, CMS-CSST-2021-A06).WNB acknowledges support from NSF grant AST-2106990.PBH is supported by NSERC grant 2017-05983.Based on observations with the NASA/ESA Hubble Space Telescope obtained from the Data Archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.Support for Program number HST-GO-15849 was provided through a grant from the STScI under NASA contract NAS5-26555.

Table 1 .
Target Properties

Table 3 .
Final photometry, Color, Luminosity, and Stellar mass (Oke & Gunn 1982)re reported in AB magnitudes(Oke & Gunn 1982), and color refers to either B − I or B − R. The last column lists the stellar masses estimated with CIGALE to compare with our fiducial stellar masses.
are consistently calibrated to the BH−host relations in quiescent local galaxies in Kor-

Table 4 .
Best-fit Parameters of the Scaling Relations