Kinematics of the Broad-line Region of 3C 273 from a Ten-year Reverberation Mapping Campaign

Despite many decades of study, the kinematics of the broad-line region of 3C~273 are still poorly understood. We report a new, high signal-to-noise, reverberation mapping campaign carried out from November 2008 to March 2018 that allows the determination of time lags between emission lines and the variable continuum with high precision. The time lag of variations in H$\beta$ relative to those of the 5100 Angstrom continuum is $146.8_{-12.1}^{+8.3}$ days in the rest frame, which agrees very well with the Paschen-$\alpha$ region measured by the GRAVITY at The Very Large Telescope Interferometer. The time lag of the H$\gamma$ emission line is found to be nearly the same as for H$\beta$. The lag of the Fe II emission is $322.0_{-57.9}^{+55.5}$ days, longer by a factor of $\sim$2 than that of the Balmer lines. The velocity-resolved lag measurements of the H$\beta$ line show a complex structure which can be possibly explained by a rotation-dominated disk with some inflowing radial velocity in the H$\beta$-emitting region. Taking the virial factor of $f_{\rm BLR} = 1.3$, we derive a BH mass of $M_{\bullet} = 4.1_{-0.4}^{+0.3} \times 10^8 M_{\odot}$ and an accretion rate of $9.3\,L_{\rm Edd}\,c^{-2}$ from the H$\beta$ line. The decomposition of its $HST$ images yields a host stellar mass of $M_* = 10^{11.3 \pm 0.7} M_\odot$, and a ratio of $M_{\bullet}/M_*\approx 2.0\times 10^{-3}$ in agreement with the Magorrian relation. In the near future, it is expected to compare the geometrically-thick BLR discovered by the GRAVITY in 3C 273 with its spatially-resolved torus in order to understand the potential connection between the BLR and the torus.


INTRODUCTION
After its discovery (Schmidt 1963), 3C 273 has been intensively studied over the whole range of the electromagnetic spectrum (from radio to γ-ray bands) and has become one of the most representative quasars (Courvoisier 1998), as it shows most of the intriguing features of active galactic nuclei (AGNs), e.g., the jet and the large flux variations at all wavelengths. Although 3C 273 is classified as a blazar for its high-energy radiation above 100 MeV, it shows a prominent big blue bump and strong emission lines in the UV and optical bands. These coexistent characteristics make 3C 273 an interesting object (Türler et al. 1999). However, the ge-ometry and kinematics of its broad-line region (BLR) are far from fully understood.
It has been shown that reverberation mapping (RM; see, e.g., Bahcall et al. 1972;Blandford & McKee 1982;Peterson et al. 1993Peterson et al. , 1998Peterson et al. , 2002Peterson et al. , 2004Kaspi et al. 2000Kaspi et al. , 2007Bentz et al. 2008Bentz et al. , 2009Denney et al. 2009;Barth et al. 2011Barth et al. , 2013Barth et al. , 2015Rafter et al. 2013;Du et al. 2014Du et al. , 2015Du et al. , 2016bDu et al. , 2018bWang et al. 2014;Shen et al. 2016;Jiang et al. 2016;Fausnaugh 2017;Grier et al. 2017;De Rosa et al. 2018) is a reliable way of measuring the masses of black holes (BHs) in AGNs (but time-consuming for very massive black holes) and investigating the geometry and kinematics of their BLRs. In Kaspi et al. (2000), the light curves of the Hα, Hβ, and Hγ emission lines of 3C 273 versus the variation of its continuum at 5100Å yielded slightly different time lags, and the cross-correlation functions (CCFs) peak at lags of τ Hα ∼ 500 days, τ Hβ ∼ 380 days, and τ Hγ ∼ 300 days, respectively. The other emission lines, such as Fe II, have not been reported yet. The large season gaps and the low sampling cadence (about 30 nights) of the light curves of 3C 273 in Kaspi et al. (2000) plausibly influence the accuracy of its time lag measurement, and result in fairly large uncertainties. Moreover, the velocity-resolved RM, which measures the time lags of the emission line at different velocities, is gradually adopted to understand the BLR kinematics and structure in recent years, and has been applied to more than a dozen AGNs (e.g., Bentz et al. 2008Bentz et al. , 2009Bentz et al. , 2010Denney et al. 2009Denney et al. , 2010Grier et al. 2013;Kollatschny & Zetzl 2013;De Rosa et al. 2015;Du et al. 2016a;Lu et al. 2016;Pei et al. 2017;Xiao et al. 2018). But higher quality spectra are needed for velocity-resolved analysis than for just getting a mean time lag of an emission line. With the very massive black hole (10 8 ∼ 10 9 M ) in the center of 3C 273, it varies on a time scale of years (see Figure 3 in Kaspi et al. 2000). Thus, a long-term RM campaign with high calibration precision, which can continue for several years, is required.
In this paper, we report a new RM campaign of 3C 273 with high quality which spans about ten years (from Nov 2008to Mar 2018. In Section 2, we describe the observation and data reduction. The light curve measurement and the inter-calibration of the data from the different telescopes are provided in Section 3. The time lags of the emission lines and the velocity-resolved RM result are given in Section 4. The BH mass is estimated in Section 6.2, and the host decomposition of 3C 273 is also provided in this section. We briefly summarize this paper in Section 7. We adopt the ΛCDM cosmology with Ω M = 0.32, Ω Λ = 0.68, and H 0 = 67 km s −1 Mpc −1 (Planck Collaboration et al. 2018) in this paper, and the luminosity distance is 782.7 Mpc.

Steward Observatory spectropolarimetric monitoring project
The spectroscopic and photometric data of 3C 273 used in this paper are mainly derived from the Steward Observatory (SO) spectropolarimetric monitoring project 1 , which is a long-term optical program to support the Fermi γ-ray Space Telescope. The project began just after the launch of Fermi in 2008 and provides almost a decade of spectropolarimetric, photometric and spectroscopic data for more than 70 blazars (Smith et al. 2009).
The SO campaign utilizes both the 2.3 m Bok Telescope on Kitt Peak and the 1.54 m Kuiper Telescope on Mt. Bigelow in Arizona. All of the observations are carried out using the SPOL spectropolarimeter (Schmidt et al. 1992), which is 1 http://james.as.arizona.edu/ ∼ psmith/Fermi/ a versatile, high-throughput, and low-resolution spectropolarimeter. A 600 mm −1 grating is used, providing a wavelength coverage of 4000-7550Å and a spectral resolution of 15∼25Å (FWHM of the line spread function) depending on the slit width chosen (Smith et al. 2009). Four slits were used for the the observations of 3C 273: 4.1 , 5.1 , 7.6 , and 12.7 . Primarily, the 7.6 -wide slit was used for the spectropolarimetric observations, which yields high-S/N spectra of the object. The 12.7 -wide slit was used for shorter observations of the quasar and a calibrated field comparison star to determine the apparent magnitude of 3C 273 in a synthetic Johnson V filter bandpass and thereby calibrate the spectrophotometry.
There are a total of 374 spectroscopic epochs for 3C 273 up until 2018 March from the SO program. Flux calibration of the spectra is based on the average sensitivity function derived from multiple observations of spectrophotometric standard stars (BD+28 4211 and/or G191 B2B) obtained during each observing campaign (typically about a week in length). Further night-to-night relative calibration was accomplished by convolving the spectra with a standard Johnson V filter bandpass and then scaling the spectrum to agree with the Vband photometric light curve.
3C 273 has 297 V -band photometric epochs, and their magnitudes were calibrated by differential photometry using the star C in Figure 1 (see more details in Smith et al. 1985). The V -band light curve from SO is shown in the following Section 3.3.
For the 374 spectroscopic observations, there are only 291 spectra that have the photometric observations in the same nights. They were calibrated using the corresponding photometric magnitudes, and are used in the following analysis. After combining the exposures in the same night, there are totally 283 spectroscopic epochs. The other 83 spectra with no simultaneous photometric observations are not used in this paper.
It should be noted that the aperture used in the spectroscopic observation includes almost all of the starlight from the host galaxy and the host contribution in the aperture is less than 6% (estimated from the image fitting in Section 6.2.3). Thus, the night-to-night calibration should not be significantly influenced by the host contribution in the different apertures used in the spectroscopy and photometry. The consistency of the SO and Lijiang light curves in the following Section 3.3 also supports that.
2.2. SEAMBH Campaigns of the Lijiang 2.4m telescope As a candidate for super-Eddington accreting massive black holes (SEAMBHs), 3C 273 was selected as a target in our SEAMBH campaign (Du et al. , 2016b(Du et al. , 2018bWang et al. 2014) for cosmology . We monitored 3C 273 from December 28, 2016 to June 13, NOTE-F 5100 is the continuum flux at 5100Å in units of 10 −15 erg s −1 cm −2Å−1 . The emission-line flux is given in units of 10 −13 erg s −1 cm −2 . Obs. is the observatory, "S" = SO, "L" = Lijiang, and "A" = ASAS-SN. This table is available in its entirety in a machine-readable form in the on-line journal. A portion is shown here for guidance regarding its form and content.  = 11.87;Smith et al. 1985) is used as the spectrophotometric standard star for the SO observation. Star G was used as the spectroscopic comparison star for the observations from Lijiang Observatory, while stars in circle were used to calibrate the Lijiang photometry. Here, we just adopt the same historical designation of "C" and "G" in Smith et al. (1985). 2017 using the 2.4 m telescope at the Lijiang Station of the Yunnan Observatories, Chinese Academy of Sciences. It is mounted with the Yunnan Faint Object Spectrograph and Camera (YFOSC), which is a multi-functional instrument both for photometry and spectroscopy. The YFOSC detector is a 2k × 4k back-illuminated CCD with pixel size of 13.5 µm, providing a 10 × 10 field of view (0 .283 per pixel). Grism 14 and a 2 .5-width slit were used, which provide a spectral dispersion of 1.74Å pixel −1 and a resolution of ∼ 500 km/s (Du et al. 2016a). The wavelength coverage is 3800-7200Å, and the wavelength was calibrated by neon and helium lamps. The spectra were reduced using the standard IRAF v2.16 procedures.
We adopted the calibration approach based on in-slit comparison star (e.g., Maoz et al. 1990;Kaspi et al. 2000;Du et al. , 2018b rather than the [O III]-based approach (e.g., van Groningen & Wanders 1992;Peterson et al. 2013;Barth et al. 2015;Fausnaugh 2017) for the spectroscopic observation. We observed 3C 273 and a nearby comparison star (star G in Figure 1) simultaneously by rotating the long slit, and used the comparison star as a standard to do the calibration. The fiducial spectrum of the comparison star is generated by averaging the spectra observed in photometric conditions. This approach ensures highly accurate relative flux calibration with a precision of ∼ 2% (Du et al. , 2018bLu et al. 2016).
We also performed photometric observations using the Johnson V filter. The light curve of 3C 273 in V -band is measured by differential photometry using 4 stars (the circled stars in Figure 1) in the same field. There are 28 photometric and 27 spectroscopic epochs in Lijiang observations, respectively. The average sampling interval is ∼ 6 days. The photometric light curve from Lijiang, which has been scaled to match the SO V -band light curve, is shown in Section 3.3.

All-Sky Automated Survey for Supernovae project
For comparison, we also show the V -band light curve of 3C 273 from the All-Sky Automated Survey for Supernovae (ASAS-SN) project 2 . The ASAS-SN project (Shappee et al. 2014;Kochanek et al. 2017) is working toward imaging the entire sky every night down to V ∼ 17 magnitude, and spans ∼ 2 − 5 years with ∼ 100 − 400 observational epochs for different objects. The ASAS-SN light curve of 3C 273 used here consists of 348 epochs in a 6-year period after averaging multiple observations obtained within a night. The ASAS-SN V -band light curve, after scaling to the SO light curve, is provided in the following Section 3.3.

Mean and RMS spectra
Small wavelength shifts in the spectra, which average to about 1.4Å for the SO data and 1.2Å for the Lijiang data, caused by instrument flexure, uncertainties in the wavelength calibration, and any mis-centering of the object within the slit have been corrected by aligning all spectra to [O III]λ5007. In addition, correction for the Galactic extinction of A V = 0.057 has been applied to the spectra (Cardelli et al. 1989;Schlafly & Finkbeiner 2011).
To show the general properties of the spectra and evaluate the variations at different velocities, we plot the mean and the root-mean-square (RMS) spectra in Figure 2. The definitions of the mean and RMS spectra arē and respectively, with N being the total number of spectra, and F i λ being the i-th spectrum. The strong signatures of emission lines in the RMS spectrum indicate that their variations are significant.

Light curves
We use the multi-component spectral fitting approach to measure the light curves of the emission lines (e.g., Barth   al. 2013; Hu et al. 2015Hu et al. , 2016. The continuum is modeled with a power-law and the narrow emission-line components ([O III]λλ4959, 5007 and narrow Hβ) by a single Gaussian. The emission features of Fe II are modeled using the template from Boroson & Green (1992). The profile of broad Hβ line is modeled by 3 Gaussians. In addition, there are some other narrow emission lines in the optical Fe II region ( Vanden Berk et al. 2001). The high-S/N spectra of 3C 273 allow us to add extra line to improve the fitting, which can also compensate the potential inconsistency between the Fe II of 3C 273 and the Fe II template of I Zw 1 from Boroson & Green (1992). Similar to Hu et al. (2015), we add a coronal line [Fe VI]λ5176 in our fitting. All of the narrow lines are fixed to have the same width and shift. The [O III]λ4959 and the narrow Hβ are fixed to have 1/3 (Osterbrock & Ferland 2006) and 1/10 (Kewley et al. 2006;Stern & Laor 2013) of the [O III]λ5007 flux, respectively. We tried to fit the broad Hβ using two Gaussians, but got poor fitting results with significant broad Hβ signal in the residual spectra, which means two Gaussians are not enough to describe the Hβ profile. Therefore, we changed to using three Gaussians in the fitting. The contribution of the host galaxy in the slit is estimated to be 6% (estimated from our image fitting in Section 6.2.3), thus is ignored in the fitting. The fitting is performed mainly at 4430-5550Å (the Hβ and Fe II region), and the window 4170-4260Å is also added to constrain the continuum slope (Hu et al. 2008). We do not add the Hγ region in the fitting window in order to reduce the number of model components. An example of a fit is shown in Figure 3.
After subtracting the continuum, Fe II, [O III]λλ4959,5007, and [Fe VI]λ5176, the light curves of the Hβ and Hγ emis-sion lines are measured by integrating the fluxes of the residual spectra in the windows 4800-4920Å and 4280-4420Å in the rest frame, respectively. We do not subtract narrow Hβ in the integration window; thus, the Hβ light curve contains the contribution from its narrow component. Considering that the variability time scale of narrow line is much longer, it does not influence the lag measurement. The 5100Å continuum and the Fe II (4434-4684Å) light curves are derived directly from the fitting results.
Spectra obtained in very poor weather conditions, or at extremely large airmasses are removed from the analysis. This includes 11 SO spectra and three spectra from Lijiang. In addition, some data points in the V -band light curve were edited out of the analysis for similar reasons and are colored gray in Figure 4. Error bars shown in the light curves contain both Poisson noise and systematic uncertainties 3 . The systematic uncertainties are estimated by the median filter method . We first smooth the light curve by a median filter with 5 points, and then subtract the smoothed light curve from the original one. The standard deviation of the residual is used as the estimate of the systematic uncertainty.
3C 273 emits strong radio and high-energy radiation (e.g., Courvoisier 1998;Türler et al. 1999;Soldi et al. 2008). Fortunately, the contribution from the synchrotron light in the optical band is weak (Impey et al. 1989;Smith et al. 1993), which is also supported by the non-blazar-like flux variations in its optical band. Meanwhile, the SO polarimetry shows that its polarization (potentially produced by the synchrotron emission) remains exceptionally low, except at the beginning of the program.
Furthermore, the RMS spectrum in Figure 2 shows stronger signal in the blue part, which means 3C 273 is more variable at shorter wavelengths. This is also evidence for the dominant contribution from the accretion disk in the optical band, because the non-thermal emission from the jet is much redder. The photometric monitoring campaign from Zeng et al. (2018) also confirms this result. In addition, the variability analysis in Chidiac et al. (2016Chidiac et al. ( , 2017 did not show any correlation between the optical/IR and the radio/X-ray/γ-ray emission, which also suggests the optical flux is dominated by the thermal emission from the accretion disk. Only a few flare-like events are found in the V -band light curve (e.g., Julian date 100 and 580, from the zero point of 2454700 in Figure 4), but they do not cause serious disruption to the time-series analysis. Therefore, the jet emission of 3C 273 does not influence the lag measurements in the following sections.

Inter-calibration
Because of the different apertures used for the SO, Lijiang, and ASAS-SN observations, inter-calibration of the three light curves is necessary. We adopt the method used by Peterson et al. (1991) to scale the Lijiang light curve to match the SO light curve. For the continuum, we assume where G 0 is a constant to account for the difference in the contribution of the host galaxy (aperture effect), ϕ 5100 is a scale factor, F S (5100Å) and F L (5100Å) are the 5100Å continuum fluxes of adjacent epochs (< 4 day apart) of the SO and Lijiang observations. For the Balmer lines and Fe II emission, we assume where "line" refers to the relevant emission line, ϕ line is the scale factor, and F S and F L are the emission-line fluxes measured for epochs < 4 days apart. Considering that the central wavelength of Hβ is not far from 5100Å, we assume ϕ Hβ = ϕ 5100 . Using the Levenberg-Marquardt algorithm (Press et al. 1992), the best-fit values are ϕ 5100,Hβ = 0.9986, ϕ Hγ = 1.0571, ϕ FeII = 0.9906, and G 0 = 0.02 × 10 −15 erg s −1 cm −2Å −1 . The small value of G 0 indicates that the aperture effect with respect to the host galaxy contribution to the spectrum is negligible and that the overall contamination by the host galaxy to the data is weak (see Section 2.1). The calibrated light curves are shown in Figure 4 and provided in Table 1. Similarly, the ASAS-SN and Lijiang V -band light curves shown in Figure 4 have also been scaled to match the SO V -band light curve.

Variability characteristics
The F var parameter (Rodríguez-Pascual et al. 1997;Edelson et al. 2002) is adopted to quantify the variability amplitude of the light curve. It is defined as where f is the mean flux in the light curve, σ is the square root of the variance and ∆ 2 is the mean square value of the uncertainties ∆ i  Here Fe II is the blue blend from 4430Å to 4680Å. The gray dashed line in the 5100Å light curve is the linear fit for detrending. The black points in panels b-f are obtained from Lijiang after the inter-calibration (see Section 3.3). All of the other spectroscopic data points are from the SO observations. In the right panels, black curves denote the CCFs (or ACFs). Blue histograms show the corresponding CCCD. The CCFs are obtained from F 5100 (after detrending) and the corresponding emission-line light curves. The January 1 of each year is marked on the time axis of the upper panel.
The F var of each light curve is listed in Table 2. The variation amplitudes of the continuum and Hβ fluxes are 10.9% and 5.2%, respectively. In addition, we also list R max in Table 2, which is defined as the ratio of the maximum to the minimum in the light curve, as well as the median values of the light curves and their standard deviations.

Time lags
The continuum light curve of 3C 273 shows long-term variations that do not appear in the emission-line light curves (see Figure 4). This long-term trend may bias the lag detection, and thus we "detrend" the 5100Å light curve by subtracting a linear fit before the lag measurement (suggested by Welsh 1999, see also, e.g., Peterson et al. 2004;Denney et al. 2010). This fit is shown in Figure 4. The detrending improves the correlation coefficients in the following crosscorrelation analysis. We tried to detrend the emission-line light curves, but the slopes of their linear fits are all close to NOTE-The emission-line fluxes and 5100Å continuum are in units of 10 −13 erg s −1 cm −2 and 10 −15 erg s −1 cm −2Å−1 , respectively. The flux, Fvar, and Rmax of the F 5100 light curve are measured in the observed frame (with Galactic extinction correction).
zero within 3-σ uncertainties. Therefore, we only detrend the 5100Å light curve. We adopt the interpolated cross-correlation function (ICCF; see Gaskell & Sparke 1986;Gaskell & Peterson 1987;White & Peterson 1994) to measure the time lags of the emission lines relative to the variation of the continuum at 5100Å after detrending. The centroid of the CCF above a typical value (80%) of the maximum correlation coefficient (r max ) is used as the lag measurement.
The uncertainties of the time lags are obtained through the "flux randomization/random subset sampling (FR/RSS)" method (Maoz et al. 1990;Peterson et al. 1998Peterson et al. , 2004, which both randomizes the measured fluxes based on their uncertainties and resamples the light curves. The resampling is performed by randomly selecting N points from the light curve (redundant selection is allowed), where N is the total number of epochs in the light curve. A subsample is constructed after ignoring the redundant points (roughly a fraction of ∼ 0.37). We repeat this process 10000 times and obtain the cross-correlation centroid distribution (CCCD) by performing CCF to the light-curve subsamples. It is similar to the so-called "bootstrapping" technique. Then the values of the time lags and the corresponding uncertainties are obtained from the median, 15.87%, and 84.13% (1σ) quantiles of the CCCD.
The auto-correlation functions (ACFs), CCFs, and CCCDs between the emission-line light curves and the detrended 5100Å light curve are shown in the right panels of Figure  4. We also provide the results without the detrending in Appendix. The time lags of the emission lines are listed in Table  3. The time lag of the Hβ emission line is 146.8 +8.3 −12.1 days in the rest frame, consistent with the lag found for the Hγ line. The response of Fe II to continuum variations is found to be nearly a factor of 2 longer than that of the Balmer lines, at 322.0 +55.5 −57.9 days in the rest frame. The lag uncertainties in our measurements are smaller than those in Kaspi et al. (2000), because of our denser sampling cadence, higher S/N ratios, and longer monitoring period. Hu et al. (2015) show a tentative correlation between the time lag ratio (τ Fe /τ Hβ ) and the intensity ratio (F Fe /F Hβ ) of the Fe II and Hβ lines. The Fe II lag is similar to the Hβ lag if F Fe /F Hβ 1, while the Fe II lag tends to be longer if F Fe /F Hβ < 1. The Fe II and Hβ time lags of 3C 273 are in good agreement with the τ Fe /τ Hβ -F Fe /F Hβ correlation ( Figure 5).

Velocity-resolved lags of Hβ
We use velocity-resolved RM to investigate the geometry and kinematics of the BLR in 3C 273. The Hβ line is divided into 15 velocity bins, with each bin having the same flux as in the RMS spectrum created after the subtraction of continuum, Fe II, and narrow lines in Section 3.2 (e.g., Bentz et al. 2008;Denney et al. 2010;Grier et al. 2013;Du et al. 2016a;De Rosa et al. 2018). The fluxes in each bin are integrated to generate light curves that are used to cross-correlate with the 5100Å variations (after detrending). The time lag and the associated uncertainties of each bin are obtained using the same method described in Section 4.2.
We plot the velocity-resolved time lags (in the observed frame) in Figure 6. The maximum correlation coefficients for individual bins are also shown in Figure 6. The velocityresolved result shows a complex structure. In general, the lags are longer at small velocities and shorter at high velocities, which is the signature of a Keplerian disk or virialized motion. However, the lags at the blue velocities are  systematically longer than those in the red wing. There are at least several possibilities: (1) a classic explanation is that the Hβ-emitting gas has some inflowing radial velocity (e.g., Bentz et al. 2008;Denney et al. 2010;Grier et al. 2013;Du et al. 2016a;De Rosa et al. 2018); (2) a rotating disk wind model suggested by Mangham et al. (2017) can produce the velocity-resolved lags resembling the "red-leadsblue" inflow signature; and (3) the BLR gas is located in an eccentric or lopsided disk, and shows net velocity along the line of sight. The BLR dynamical modeling technique (e.g., Pancoast et al. 2011Pancoast et al. , 2012Grier et al. 2017;Williams et al. 2018) is required to investigate the geometry and kinematics of the BLR of 3C 273 in more detail in the future.

Line widths
We investigate the behavior of the broad emission-line profiles of Hβ and Hγ through the measurement of the FWHM and the line dispersion (σ line ). The line dispersion is defined as where λ 0 = λf (λ)dλ/ f (λ)dλ is the flux-weighted centroid wavelength of the emission line profile, λ 2 = λ 2 f (λ)dλ/ f (λ)dλ, and f (λ) is the flux density. Line widths are measured from the line profiles in the mean and RMS spectra after subtraction of the continuum, Fe II emission, and narrow emission lines. The uncertainties of the line widths are estimated by the bootstrap method (Peterson et al. 2004). We randomly select N spectra with replacement (redundant selection is allowed, N is the total number of spectra for 3C 273), and construct a spectral subset after removing the redundant spectra. Then, the mean and RMS spectra are generated from the subset, and are used to measure the line widths. This process is repeated 10000 times to create the line width distributions. The standard deviations of the FWHM and σ line distributions are regarded as the uncertainties.
The line widths of Hβ and Hγ are calculated separately for the SO and Lijiang data. To estimate the spectral broadening caused by the instruments and seeing, we compare our spectra with higher-resolution spectra obtained with an instrument with a well-understood line spread function. In this case, we use the o44301010 and o44301020 observations of 3C 273 made by STIS aboard the Hubble Space Telescope (Hutchings et al. 1998). The spectral broadening of SO and Lijiang are 1157 ± 46 km/s and 485 ± 38 km/s, respectively. After subtracting the spectral broadening in quadrature, the line widths of the Hβ and Hγ lines are listed in Table 4. For the SO spectra, the FWHM of Hβ and Hγ are both around ∼3300 km/s. We adopt the FWHM from the mean spectrum to calculate of BH mass in the following sections, because the σ line measurement is sensitive to the wings of the emission line and the accuracy of the continuum subtraction. The FWHM of Fe II is 2142.4±116.1 km/s, which is obtained directly from the fitting. Compared with the Hβ and Hγ, the FWHM of Fe II is significantly narrower. This suggests that the Fe II-emitting region is farther away from the BH than the region in the BLR producing the Balmer emission.
The FWHM and σ line of individual spectra are also measured to investigate the line width changes during the campaign. After correcting for spectral broadening, the FWHM and σ line of Hβ from SO for the entire monitoring period are shown in Figure 7. Figure    σ line vs. F Hβ , respectively.) To quantify the correlation, we perform the linear regression of where V Hβ is FWHM Hβ or σ Hβ . The uncertainties of the parameters are obtained using the bootstrap technique. The results are shown in Figure 8. The best-fit values for (α, β) are (−1.28 ± 0.37, −0.41 ± 0.03) and (1.46 ± 0.15, −0.15 ± 0.01) for flux versus FWHM and σ Hβ , respectively. Similar correlations are reported for some other objects (e.g., Barth et al. 2015).

Line width versus lags
The RM campaign of 3C 273 performed by the International Ultraviolet Explorer (IUE) yielded time lags of τ Lyα+Nv = 435 +77 −77 and τ C IV = 690 +102 −106 days in the rest frame for Lyα + N V and C IV lines relative to the UV continuum variation (Paltani & Türler 2005). The echo of the Lyα + N V lines with respect to the UV continuum is quite good with r max ≈ 0.6, but the C IV response remains largely uncertain. Paltani & Türler (2005) took the centroid lag above 0.6 r max of the Lyα + N V CCF as the lag measurement. For comparison to our optical measurements, we re-calculated their Lyα + N V lag using the 0.8 r max criterion. There are multiple peaks in the ultraviolet CCF, and we only adopt the highest peak in the calculation. The re-determined Lyα + N V lag is τ Lyα+Nv = 72.2 +58.1 −23.5 days in the rest frame. We plot the correlation between the lags and the FWHMs both for the UV and optical lines in Figure 9 , where it is seen that FWHM ∝ τ −0.55±0.01 for the optical emission lines. This is similar to the results for NGC 5548 (Peterson & Wandel  1999). This relation steepens to FWHM ∝ τ −1.1±0.3 with the inclusion of the Lyα + N V lines. 3C 273 shows a FWHM ∝ τ −0.56 relation for its optical emission lines that the UV emission lines do not follow. It implies that the geometry or kinematics of the UV-emitting clouds are different from the optical-line region, i.e., their virial factors in the BH mass measurements should be very different. The optical emission lines are relatively closer to FWHM ∝ τ −0.5 relation (consistent with the velocityresolved result that the kinematics in the Hβ region is generally a Keplerian disk or virialized motion); thus, we prefer to estimate the BH mass based on those lines. Given that the IUE campaign was carried out more than a decade before the optical RM program, the steeper slope may also be caused by a physical change in the geometry and/or kinematics of the BLR.

BLACK HOLE MASS AND ACCRETION RATE
Combining the results for the time lag, τ line , and the line width, V line , the BH mass can be estimated by where R BLR = cτ line is the emissivity-weighted radius of the BLR, V line is given by the FWHM or σ line , c is the speed of light, G is the gravitational constant, and f BLR is the virial factor depending on the kinematics, geometry, and inclination angle of the BLR. In general, the average f BLR can be estimated by comparing RM-based BH masses for AGNs with available bulge stellar velocity dispersions with BH masses predicted from the M − σ relation of inactive galaxies (e.g., Onken et al. 2004;Ho & Kim 2014;Woo et al. 2015;Park et al. 2012;Grier et al. 2013;Ho & Kim 2014). Specifically, for FWHM of Hβ measured from their mean spectra, Ho & Kim (2014) demonstrate that the f BLR is 1.3 for AGNs with classical bulges or elliptical host galaxies (see also Mejía-Restrepo et al. 2018). Independent from the calibration based on the M • − σ * relation, the Bayesian-based BLR dynamical modeling method (Pancoast et al. , 2012 has been adopted in recent years to provide f BLR for individual AGNs with different BLR geometry and kinematics (e.g., Pancoast et al. 2014aPancoast et al. ,b, 2018Grier et al. 2017;Li et al. 2013Li et al. , 2018Williams et al. 2018). The velocity-resolved time lags of 3C 273 show "red-leadsblue" signatures, which may arise from inflow kinematics (in Section 4.3). Virial factors for AGNs with inflowing signatures as derived from dynamical modeling are diverse and range from 0.58 − 3.98 (Grier et al. 2017;Pancoast et al. 2014b), but they average close to the value found by Ho & Kim (2014). Thus, we simply adopt f BLR = 1.3 to calculate the BH mass for 3C 273, but acknowledge that there is a large amount of uncertainty in the value of this parameter (we do not include the uncertainty of f BLR in our mass calculation). Based on the measurements of τ Hβ = 146.8 +8.3 −12.1 days (in the rest frame) and FWHM Hβ = 3314.1 ± 59.3 km/s, the BH mass of 3C 273 is M • = 4.1 +0.3 −0.4 × 10 8 M . The relativistic jet is found to have an inclination of 10 • with a Lorentz factor of Γ ≈ 11 (Davis et al. 1991;Abraham & Romero 1999), which is nicely agreeable with the measurements (i = 12 • ± 2 • ) of the GRAVITY (see Section 6.1), showing that the BLR and the jet are perpendicular to each other (Gravity Collaboration et al. 2018). Since Θ BLR ≈ 45 • , the inclination (12 • ) implies that the present estimate of SMBH mass is less affected by inclinations since the normalization factor f BLR ∝ (H BLR /R) 2 + sin 2 i −1 is dominated by the geometric factor of H BLR /R (i.e., H BLR /R = tan Θ BLR , where H BLR is the half-thickness of the BLR). The GRAVITY observations provide an SMBH mass of M • = (2.6 ± 1.1) × 10 8 M (Gravity Collaboration et al. 2018), which is consistent with the present results within the uncertainties. We would like to point out that the most uncertainties of SMBH masses from Eq. (10) are due to f BLR . Detailed modelling BLR through Markov Chain Monte Carlo (MCMC) simulations can be done for the mass as well as the BLR geometry (Pancoast et al. , 2014aLi et al. 2013Li et al. , 2018 independent of this factor. Application of the MCMC NOTE-"Obs. season" cadence means that it is averaged only during observable seasons and the entire one does that it is averaged for the entire campaign. Length refers to the total period from the start to end. Nspec is the number of spectra taken from the corresponding campaigns. technique 4 to 3C 273 will be carried out separately (Li et al. 2019 in preparation). We use the equation for the dimensionless accretion rate oḟ the value corresponding to the cosmological parameters used in this paper), we find that L bol /L Edd ≈ 2.4, indicating that the BH is accreting in this quasar at a mildly super-Eddington rate.

Comparing with the previous and the GRAVITY results
Before the present work, Kaspi et al. (2000) made a mapping measurement of 3C 273 through a 7-year campaign of the Bok telescope and the Wise 1m telescope. They obtained the centroid lags of reverberation of the Balmer lines, (R Hα , R Hβ , R Hγ ) = (514 +65 −64 , 382 +117 −96 , 307 +57 −86 )ltd, respectively, that are twice more as long as the present results. But the error bars of the present measurements are much smaller than those given by Kaspi et al. (2000). Peterson et al. (2004) revisited the data of Kaspi et al. (2000) and found R Hβ = 306.8 +68.5 −90.9 ltd, it is still double the present measured result in this paper. A brief comparison of the two campaigns is given by Table 5, showing that the present campaign has a much higher cadence than Kaspi et al. (2000). Considering the BLR dynamical timescales of τ BLR ≈ R BLR /V FWHM = 37.4 R 150 V −1 3300 yrs, the BLR cannot change by a factor of 2 in the last 27 yrs (since 1991), where R 150 = R BLR /150ltd and V 3300 = V FWHM /3300km s −1 . Moreover, variations of the 5100Å luminosity are about ∆L/L 5100 1/4 during the two campaigns and thus the difference of Hβ lags measured by the campaigns is not caused by luminosity changes. The present measurements are more robust and reliable.
Very recently 5 , the GRAVITY instrument installed on the Very Large Telescope Interferometer (VLTI) found a flatten disk structure of the BLR in 3C 273 with a half-opening angle of Θ BLR = 45 •+9 −6 and an emissivity-averaged radius of the Paschen-α emission region R Paα = 145 ± 35ltd (Gravity Collaboration et al. 2018). It is exciting to find that R Paα = R Hβ = R Hγ holds within their error bars as listed in Table  3. This lends a robust demonstration that RM measurements are very reliable to spatially resolve compact regions in time domain through small telescopes.
Moreover, the opening angle of the geometrically thick BLR in 3C 273 might agree with that of the dusty torus. This is helpful to understand the hypothesis that the BLRclouds are supplied by tidally captured clumps from the torus (Wang et al. 2017). In principle, the torus opening angles can be estimated by infrared emission as the reprocessed emission of optical and UV luminosities (Wang et al. 2005;Cao 2005;Maiolino et al. 2007), but the infrared emissions in 3C 273 could be seriously contaminated by non-thermal emission from the jet. So it's hard to estimate the opening angle of the torus in this way. Fortunately, the [O III] luminosity is L [O III] ≈ 1.5 × 10 43 erg s −1 from this campaign, allowing an estimate of the torus half-opening angle of Θ torus ≈ 30 • ± 10 • from the correlation between type 2 quasar ratio and [O III] luminosity (i.e., the receding torus model, e.g. Fig. 12 in Reyes et al. 2008). This shows potential evidence for the hypothesis of the physical connection between the BLR and the torus. The GRAVITY observation of the torus in 3C 273 in the future, as done in NGC 1068 (Jaffe et al. 2004), will directly measure the torus geometry to check if the geometrically-thick BLR matches that of the dusty torus. This will observationally justify the hypothesis of the BLR origin from the failed outflows of the outer part accretion disk (Czerny & Hryniewicz 2011) or from torus (Wang et al. 2017 AGNs show strong correlations with their host galaxies, for example, the BH mass is strongly correlated with the mass of the galactic bulge of the host (e.g., Magorrian et al. 1998; Kormendy . To derive the bulge mass in 3C 273, we analyze its HST/WFC3 images 6 , which were observed on March 17, 2013 (GO-12903, PI: Luis C. Ho) with the F547M filter in the UVIS channel for 360 sec and with F105W filter in the IR channel for 147 sec. Two additional short F547M exposures for 18 sec and 6 sec were taken to warrant against saturation of the nucleus in the long exposure. These two filters were selected to mimic the B and I bands in the rest frame of 3C 273. To better sample the point-spread function (PSF), the long UVIS observation was taken with a threepoint linear dither pattern, while the IR observation was taken with the four-point box dither pattern. To avoid overheads due to buffer dump, we employed the UVIS2-M1K1C-SUB 1k×1k subarray for the UVIS channel and the IRSUB512 subarray for the IR channel, which results in a restricted fieldof-view of 40×40 arcsec and 67×67 arcsec, respectively. Because of the sparseness of field stars in the vicinity of 3C 273, the observations were conducted in GYRO mode, leading to a considerable degradation of the PSF of the dither-combined image generated from the standard data reduction pipeline. Instead, we use the DrizzlePac task As-6 There is another archive data of HST observation led by Bentz et al. (2009), but only provided one color image. We revisited her data and obtained similar results with Bentz et al. (2009) to remove host contaminations. However, we need two-color images of HST observations for stellar mass in this paper. troDrizzle (v1.1.16) to correct the geometric distortion, align the sub-exposures, perform sky subtraction, remove cosmic rays, and, finally combine the different exposures. The core of the AGN was severely saturated in the long F547M exposure, and it was replaced with an appropriately scaled version of the 6 sec short exposure. Because the PSF of HST/WFC3 is under-sampled, we broaden the science image and the PSF by a Gaussian kernel such that the final images are Nyquistsampled (Kim et al. 2008a). This largely removes the subpixel mismatch of the core of the PSF.

Image decomposition
The F105W filter image for 3C 273 shows a dominant nucleus and a host galaxy with no clear disk component. The F547M filter image is considerably shallower. To extract quantitative measurements of the bulge, we use the program GALFIT (Peng et al. 2002(Peng et al. , 2010 to fit two-dimensional surface brightness distributions to the HST images. A crucial ingredient is the PSF, which will have a strong effect on the brightness of the active nucleus. Unfortunately, no suitable bright star is available to be used as the PSF within the limited field-of-view of the subarray WFC3 images. Instead, we generated a high-S/N empirical PSF by combining a large number (24 for F105W and 12 for F547M) of bright, isolated, unsaturated stars observed in other WFC3 programs. Extensive tests, consisting of fits to isolated bright stars, indicate that our stacked empirical PSF is far superior to syn- thetic PSFs generated from the TinyTim program (Krist & Hook 2008), and it has higher S/N than the PSFs of individual stars. The reduced chi-squared of the fits are ∼3 times larger for the TinyTim synthetic PSF. Comparison of empirical PSFs observed from different programs indicate that the WFC3 PSF does not vary significantly with time (< 10%).
We first obtain the best global fit on the deeper F105W image, whose redder wavelength coverage is more sensitive to the host. Models consisting of two components are adopted, with a point source (represented by the PSF) for the nucleus, and the bulge parameterized by a Sérsic function (Sérsic 1968) with index n. Models with n ≈ 3 − 5 give the best fitting residuals. We adopt n = 4 as the best model and include the difference between the n = 3 and n = 5 model into the uncertainty. Figure 10 summarizes the model fits in F105W and F547M filter for the n = 4 model. We use the m = 1 Fourier mode, which is sensitive to lopsidedness, to gauge the degree of global asymmetry of the galaxy (see, e.g., Kim et al. 2008bKim et al. , 2017. Since the host is considerably weaker in the F547M image, it is fit by keeping the structural parameters fixed to the values obtained from the F105W model, solving only for the magnitudes of the nucleus and host.
The best-fit parameters are given in Table 6. For bright type 1 AGNs, such as 3C 273, the uncertainties are dominated by errors in the PSF adopted for the nucleus. We estimate the uncertainties in the PSF by generating variants of the empirical PSF by combining different subsets of stars, and then repeating the fit. Subtraction of the background can also affect the fit. We subtract ±1σ of the background level from the image, and then take the mean of the differences relative to the fiducial value of the sky as the uncertainty. The uncertainties caused by the PSF, background subtraction, and the Sérsic index n are all taken into account in the final error estimation.

Stellar mass
The GALFIT decomposition of the HST images yields F547M and F105W magnitudes for the bulge of the host galaxy of 3C 273. We convert these HST-based magnitudes into intrinsic, rest frame I-band magnitude and B − I color, from which we can estimate the stellar mass following the prescriptions of Bell & de Jong (2001). To estimate the Kcorrection, we use Bruzual & Charlot (2003) models with solar metallicity, a Chabrier (2003) stellar initial mass function (IMF), and an exponentially decreasing star formation history with a star formation timescale of 0.6 Gyr to generate a series of template spectra with ages spanning 1 to 12 Gyr. After accounting for Galactic extinction and redshift, we convolve the spectra with the response functions of the HST filters to generate synthetic F547M and F105W magnitudes.
The observed host color is F547M − F105W = 1.6 ± 0.9. This color is best matched with a stellar population template of 2.5 Gyr (the templates of 1 Gyr and 12 Gry for the lower and upper limits), resulting in M I = −24.3 ± 0.6 and B − I = 1.5 ± 1.0. Following Bell & de Jong (2001) and the correction of Longhetti & Saracco (2009), we derive a stellar mass of M * = 10 11.3±0.7 M , assuming solar metallicity and a Chabrier IMF. The ratio of BH mass and spheroid is about 2.0 × 10 −3 , which agrees with the Magorrian relation (e.g., Kormendy & Ho 2013).

SUMMARY
We present the results from a reverberation mapping campaign of 3C 273 from Nov 2008 to Mar 2018. Time lags relative to the observed continuum variations of the active nucleus of several emission lines were successfully detected, and we measure the velocity-resolved lags for the Hβ line.
• The Fe II lines have a lag of τ Fe = 322.0 +55.5 −57.9 days in the rest frame, which follows the τ Fe /τ Hβ − F Fe /F Hβ correlation found in Hu et al. (2015). The Fe II and the Balmer line regions follow the virialized relation of τ ∝ V −2 FWHM , showing the stratified structures in space.
• The velocity-resolved lag measurements of the Hβ line show a complex structure. The lags are longer at small velocities and shorter at high velocities, which is the signature of a rotation-dominated disk. In addition, the blue wing shows longer lags than those at red velocities. This may be explained by inflowing radial motions in the Hβ-emitting gas or some other special BLR kinematics.
• Along with the UV lines, we find that 3C 273 has a stratified structure in its BLR, with higher-ionization lines arising from inner regions but lower-ionization lines from the outer part. The UV lines deviate significantly from the FWHM ∝ τ −0.5 relation.
• Adopting a virial factor of f BLR = 1.3, 3C 273 has a BH mass of M • = 4.1 +0.3 −0.4 × 10 8 M and is accreting with a rate of 9.3 L Edd c −2 , indicating that the black hole is undergoing intermediate super-Eddington accretion.
• From decomposition of its HST images, we obtain a host stellar mass of M * = 10 11.3±0.7 M . We find that a ratio of BH mass and host spheroid (2 × 10 −3 ) follows the Magorrian relation.
Considering that the geometrically thick BLR of 3C 273 detected by the GRAVITY is roughly consistent with the torus obtained by the receding statistic, we expect to compare with torus structure spatially resolved by the VLTI in order to reveal the origin of BLR clouds, which is suggested to be supplied by the central black hole tidal capture of clumps from the torus (Wang et al. 2017).
We thank an anonymous referee for a helpful report. We acknowledge the support by National Key R&D Program of China (grants 2016YFA0400701 and 2016YFA0400702), by NSFC through grants NSFC-11873048, -11833008, -11573026, -11473002, -11721303, -11773029, -11833008, -11690024, and by Grant No. QYZDJ-SSW-SLH007 from the Key Research Program of Frontier Sciences, CAS, by the Strategic Priority Research Program of the Chinese Academy of Sciences grant No.XDB23010400. Data from the Steward Observatory spectropolarimetric monitoring project were used. This program is supported by Fermi Guest Investigator grants NNX08AW56G, NNX09AU10G, NNX12AO93G, and NNX15AU81G. We also acknowledge the support of the staff of the Lijiang 2.4m telescope. Funding for the telescope has been provided by CAS and the People's Government of Yunnan Province.

APPENDIX
For completeness, we also show the results of the time-series analysis and the velocity-resolved lag measurement without detrending in Figures 11 and 12. The time lags of the Hβ, Hγ, and Fe II lines are listed in Table 7. The long-term variations in 5100Å light curve (i.e., the long-term trending) bias the lag measurements and the velocity-resolved results significantly. Comparing with Table 3, we find that CCF correlations are also stronger than the results without detrending (see Table 7). We adopt the detrended results given in Table 3 of the main text.