Reverberation Mapping of Two Luminous Quasars: The Broad-line Region Structure and Black Hole Mass

We spectroscopically monitored two quasars, PG 0923+201 and PG 1001+291, both with an optical luminosity larger than 10⁴⁵ erg s⁻¹. Their basic properties are summarized below.

PG 0923+201 has a redshift of z = 0.1921, m_g = 15.6 mag, and V-band extinction A_V = 0.116 mag (Schlafly & Finkbeiner 2011). Its spectrum shows prominent Fe ii emission lines and a weak [O iii] λ5007 line (Boroson & Green 1992). The ratios of the equivalent width of [O iii] λ5007 and Fe ii between 4434 Å and 4684 Å to that of Hβ are about 0.04 and 0.72, respectively. Such strong Fe ii and weak [O iii] emission lines are common features seen in narrow-line Seyfert 1 galaxies (NLS1s), which are generally believed to be accreting at a super-Eddington rate. However, the FWHM of the Hβ line is as large as ∼7600 km s⁻¹, much broader than NLS1s (usually narrower than 2000 km s⁻¹).

PG 1001+291 has a redshift of z = 0.3272, m_g = 15.9mag, and A_V = 0.06mag (Schlafly & Finkbeiner 2011). Similar to PG 0923+201, PG 1001+291 also shows prominent Fe ii emission lines with ${{ \mathcal R }}_{\mathrm{Fe}}$ ∼ 1.17 and a weak [O iii] line (Du & Wang 2019). However, its Hβ line has a much narrower FWHM, only about 2100 km s⁻¹. PG 1001+291 was monitored between 2015 and 2017 by the Super-Eddington Accreting Massive Black Holes (SEAMBH) campaign ¹⁶ (Du et al. 2018b, who used the alternative SDSS name for PG 1001+291, namely, SDSS J100402.61+285535.3). The Hβ lag (${32.2}_{-4.2}^{+43.5}$ days) was relatively uncertain due to the short monitoring period. Nevertheless, the Hβ lag of PG 1001+291 was found to be 0.80 dex below the R_Hβ –L₅₁₀₀ relation (Du et al. 2018b).

In our observations, spectroscopy was mainly undertaken with the Lijiang 2.4 m telescope at the Yunnan Observatories of the Chinese Academy of Sciences. It is equipped with the Yunnan Faint Object Spectrograph and Camera (YFOSC), which can switch quickly between spectroscopy and photometry modes. The observations and data reduction were carried out following Du et al. (2014, 2015). We briefly summarize important points about the observation setup below. (1) For PG 0923+201, the spectra were obtained by Grism 14 with a dispersion of 1.78 Å pixel⁻¹ and a wavelength coverage of 3800–7400 Å. We used a long slit with a projected width of 25. The instrumental broadening was about 500 km s⁻¹ (in terms of FWHM; Du et al. 2014). (2) For PG 1001+291, due to its larger redshift, Grism 14 would induce contamination of the second-order spectrum at the observed-frame wavelength longer than 6600 Å (see Lu et al. 2019; Feng et al. 2021). We therefore used Grism 3 to avoid contamination. Grism 3 has a spectral coverage of 3800–9000 Å and a dispersion of 2.9 Å pixel⁻¹. We adopted a 50 wide slit and the instrumental broadening is about 1200 km s⁻¹ (Du et al. 2015).

For each target, we rotated the slit to cover a selected comparison star simultaneously for flux calibration (see Figure 1). We selected the comparison star WISEA J092607.03+195357.5 for PG 0923+201 and WISEA J100406.45+285631.6 for PG 1001+291. The spectroscopic images were reduced using standard IRAF procedures. A uniform aperture of 85 and a background region of 74–14'' were used for spectrum extraction. Two consecutive 1200 s exposures were taken on each observing night. For PG 0923+201, we obtained a total of 89 epochs of spectroscopic observations between 2017 October and 2020 May. For PG 1001+291, we obtained a total of 164 epochs between 2015 November and 2019 May. The median sampling intervals of PG 0923+201 and PG 1001+291 are 5.9 and 4.9 days, respectively (see Table 1). For each night, the typical signal-to-noise ratio (S/N) per pixel at 5100 Å (rest frame) is ∼57 for PG 0923+201 and ∼48 for PG 1001+291.

Zoom In Zoom Out Reset image size

Since 2017 May, PG 1001+291 was also observed with the Centro Astronómico Hispano-Alemán (CAHA) 2.2 m telescope at the Calar Alto Observatory in Spain (see Hu et al. 2020a, 2020b). The spectra were taken with the Calar Alto Faint Object Spectrograph (CAFOS) using the same observation strategy described above. The adopted Grism G-200 has a dispersion of 4.47 Å pixel⁻¹ and a spectral coverage of 4000–8500 Å. We used a long slit with a projected width of 30 and the resulting instrumental broadening is about 1000 km s⁻¹ (Hu et al. 2020b). We selected the same comparison star as in our Lijiang observations. A uniform aperture of 106 and a background region of 138–265 were used for spectrum extraction (see Hu et al. 2020b for the data reduction details). In total, we obtained 16 epochs between 2017 May and 2019 January at the Calar Alto Observatory. The typical S/N per pixel at 5100 Å (rest frame) is ∼68 of each night.

The Fe ii blends are relatively strong beneath the [O iii] line, and the [O iii] λ5007 line is too weak to apply [O iii]-based calibration (van Groningen & Wanders 1992). Therefore, we adopt the method based on comparison stars, which provides sufficiently accurate flux calibration (Maoz et al. 1990; Du et al. 2018b). Because the target and its comparison star are observed simultaneously, their emitted lights travel along the same path and thereby suffer the same seeing and atmospheric conditions. This enables high-accuracy relative flux measurements even in relatively poor weather conditions. The calibration steps are as follows. First, we use the spectrophotometric standard stars to calibrate the absolute spectra of the comparison star observed on several good-weather nights. We average these calibrated spectra of the comparison star to generate a fiducial spectrum. Second, in each observation, we compare the comparison star's spectrum with the fiducial spectrum to get a sensitivity function and use this sensitivity function to calibrate the target's spectrum. Below, we also illustrate that the selected comparison stars are stable in photometry and therefore suitable for flux calibration. Figure 1 plots images to illustrate sky locations of the two quasars and their comparison stars along with slits.

Photometry was obtained directly through the image mode of YFOSC. We used the Johnson V filter and the SDSS ${r}^{{\prime} }$ band filter for PG 0923+201 and PG 1001+291, respectively. Typically we took three 50 s exposures for each object on each observing night. We used standard IRAF procedures to reduce the photometric data. Several stars were selected in the $10^{\prime} \times 10^{\prime}$ field of view for differential photometry (see Figure 1). According to aperture tests, we found that a circular aperture with a radius of 45 and 42 provides the best photometry for the two objects, respectively. The inner and outer radii of the background was uniformly set to be 85 − 17''. Figure 2 shows the photometric light curves of the comparison stars over our monitoring periods, both of which have a standard deviation of ≲0.01 mag, indicating that the comparison stars were stable enough for flux calibration. We finally obtained a total of 85 and 160 epochs of photometric observations for PG 0923+201 and PG 1001+291 with median sampling intervals of 6.0 and 4.1 days, respectively (see Table 1). Tables 2 and 3 list the photometric light curves of the two targets.

Zoom In Zoom Out Reset image size

We also obtained photometry of PG 1001+291 using the CAFOS image mode with a Johnson V filter at the Calar Alto Observatory. Typically, three 60 s exposures were taken on each observing night. The photometry was reduced using standard IRAF procedures with the same configurations as in our Lijiang observations. We secured 14 epochs between 2017 June and 2019 January at the Calar Alto Observatory.

Besides our observations, we compiled photometry from archived data of the All-Sky Automated Survey for Supernovae (ASAS-SN) and the Zwicky Transient Facility (ZTF).

The ASAS-SN ¹⁷ is a long-term project designed to survey the whole visible sky every night down to about 17 mag to discover nearby supernovae and other transient sources (Shappee et al. 2014; Kochanek et al. 2017). The ASAS-SN project is composed of multiple stations, each station containing a 14 cm Nikon telephoto lenses equipped with a thermoelectrically cooled CCD camera. The field of view of each camera is about 4.5 deg², and the pixel scale is 8''. Observations generally used V- or g-band filters for three dithered 90 s exposures. The photometric data were processed using the IRAF apphot package with an aperture radius of 16'' and calibrated according to AAVSO Photometric All-Sky Survey (Henden et al. 2012). For PG 0923+201, there are 519 epochs with a median sampling interval of 1.7 days from 2015 October to 2020 June (see Table 1). For PG 1001+291, the ASAS-SN data are not adopted because of the relatively poor data quality.

The ZTF ¹⁸ is designed for transients and variables and uses the Palomar 48 inch Schmidt Telescope with a 47 deg² field of view and a 600 megapixel camera to scan the entire northern visible sky (Bellm et al. 2019; Graham et al. 2019; Masci et al. 2019). The ZTF provides two filters, ZTF-g and ZTF-r, with the median photometric depths down to 20.8 and 20.6 mag, respectively. We combine the light curves of the ZTF-g and ZTF-r bands using the intercalibration method described in Section 3.2. There are 374 epochs with typically 3 observations per night between 2018 March and 2020 May for PG 0923+201 and 143 epochs with a median sampling interval of 1.0 days from 2018 March to 2019 June for PG 1001+291 (see Table 1).

We measure the 5100 Å flux density and broad Hβ line via a multicomponent spectral fitting scheme (SFS) following Hu et al. (2015, 2020a, 2020b). Before fitting, we correct the Galactic extinction using the extinction law of Cardelli et al. (1989) with R_V = 3.1. We employ the DASpec ¹⁹ software to perform multicomponent spectral fitting. The software uses the Levenberg–Marquardt method to minimize the chi-square. The details of the fitting procedures for the two targets are slightly different, as described below.

PG 0923 + 201—There is no obvious stellar absorption feature in our spectra (see Figure 3; near ∼4940 and 5280 Å are telluric absorptions), so we do not add a host galaxy template in the fitting. The fitting components include (1) a single power law, (2) a Fe ii template from Boroson & Green (1992), (3) a single Gaussian for broad emission lines Hβ and Hγ, (4) double Gaussians for [O iii] λ λ 4959, 5007 and [O iii] λ4363 with the same tied velocity shifts and widths among the three lines. The narrow Hβ is not added because the flux of the narrow Hβ is less than 2% of the total Hβ flux, and the narrow Hβ flux is estimated by assuming a typical flux ratio of 0.1 between the narrow Hβ and [O iii] λ5007 (e.g., Veilleux & Osterbrock 1987; Hu et al. 2015).

Zoom In Zoom Out Reset image size

Because the Hγ is blended with [O iii] λ4363, we fix the velocity width and shift of Hγ to those of Hβ. The flux ratio of [O iii] λ4363 relative to [O iii] λ5007 is fixed to 0.251, obtained by the best fit of the mean spectrum. The [O iii] lines have asymmetric blue wings, which are seen more clearly in the high-resolution spectrum from the Sloan Digital Sky Survey (SDSS). We therefore apply two Gaussians to fit the [O iii] lines, one of which is set to be narrower and the other is set to be broader. The width and shift of the broader Gaussian component of [O iii] lines are fixed to the value obtained via the best fit to the SDSS spectrum after considering the different instrumental broadening of SDSS and the Lijiang telescope. The fitting is performed at 4170–5550 Å, excluding the telluric absorption windows around 4940 Å and 5280 Å.

We note that although the He ii λ4686 line is relatively prominent in the rms spectrum, it is weak and highly blended with the Fe ii and the blue wing of the Hβ in individual spectra so that it cannot be well constrained in the SFS. We therefore do not include a He ii component. To evaluate the influences of the He ii line on the time delay measurements, we present the results by adding a Gaussian to account for the He ii component in Appendix A. The obtained time delays from the two schemes are consistent within uncertainties, indicating that the influences of the He ii line are minor.

PG 1001 + 291—Again, the host galaxy template is not included. The fitting components include (1) a single power law, (2) a Fe ii template from Boroson & Green (1992), (3) a fourth-order Gauss–Hermite function for broad Hβ and broad Hγ, (4) a single Gaussian with the same velocity width and shift for narrow lines [O iii] λ λ 4959, 5007 and [O iii] λ4363. Similar to PG 0923+201, the width and shift of Hγ are fixed to those of Hβ, and the flux ratio of [O iii] λ4363 to [O iii] λ5007 is fixed to the fitted value of the mean spectrum. The fitting is performed at 4170–5390 Å, excluding the telluric absorption windows around 4730 Å and 5200 Å (see Figure 3).

The light curves of the continuum flux density at 5100 Å(F₅₁₀₀) and the broad emission lines (Hβ and Hγ of PG 0923+201; Hβ, Hγ, and Fe ii of PG 1001+291) are directly obtained from the above fitting and summarized in Tables 2 and 3. The reported errors of the light curves include the Poisson errors and additional systematic errors (added in quadrature). The additional systematic errors are determined following the procedure described in Du et al. (2014). We note that the current Fe ii data of PG 0923+201 are not enough to detect a reliable time delay so that we do not present Fe ii analysis for PG 0923+201 in this work.

We merge the photometric data from our observations, the ASAS-SN, and the ZTF into the 5100 Å continuum to obtain a combined continuum light curve. Time lags between different bands are typically of days (e.g., Edelson et al. 2019), far smaller than the time lags of broad emissions; therefore, such an intercalibration is feasible and does not affect the final time lag measurements. Due to inhomogeneous apertures, the photometric and 5100 Å continuum data need intercalibration, namely, applying additive and multiplicative factors to bring different light curves into a common scale. We use the Python package PyCALI, ²⁰ which employs a Bayesian framework to do intercalibration (see details in Li et al. 2014). The intercalibrated continuum light curves are shown in Figure 4. For PG 1001+291, the ASAS-SN data are not used (see Section 2.3). After intercalibration, we further rebin the continuum light curves by one day apart to combine measurements on the same night. The finally rebinned continuum light curves are used to measure time lags, shown in panel (a) of Figure 5 and 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For PG 1001+291, we obtained spectra from both Lijiang and CAHA observations; therefore, an intercalibration of the spectra is also required. Because the instrumental broadening widths (in terms of FWHM) of the Lijiang and CAHA observations are roughly similar and far smaller than the broad line widths, there is no need to adjust the spectral resolutions. We simply use the intercalibration factors obtained above to align the fluxes of the spectra. In panels (b)–(d) of Figure 6, we show the intercalibrated light curves of emission lines for PG 1001+291 with Lijiang data in black and CAHA data in blue.

As usual, we calculate the quantities F_var and ${R}_{\max }$ to measure the intrinsic variability amplitudes (Rodríguez-Pascual et al. 1997). Here, ${R}_{\max }$ is the ratio of the maximum to minimum fluxes of the light curve, and F_var is defined as

$$\begin{eqnarray}&&{F}_{\mathrm{var}}=\displaystyle \frac{{({S}^{2}-{\bigtriangleup }^{2})}^{1/2}}{\bar{F}},\end{eqnarray} \tag{ 2 }$$

where S² is the sample variance, △² is the mean square error, and $\bar{F}$ is the sample mean flux

$$\begin{eqnarray}&&{S}^{2}=\sum _{i=1}^{N}\displaystyle \frac{{({F}_{i}-\bar{F})}^{2}}{(N-1)},\quad {\bigtriangleup }^{2}=\sum _{i=1}^{N}\displaystyle \frac{{{\rm{\Delta }}}_{i}^{2}}{N},\quad \bar{F}=\sum _{i=1}^{N}\displaystyle \frac{{F}_{i}}{N},\end{eqnarray} \tag{ 3 }$$

where F_i is the flux of the ith observation, Δ_i is the uncertainty of F_i , and N is the total number of epochs. The standard deviation of F_var is estimated by Edelson et al. (2002):

$$\begin{eqnarray}&&{\sigma }_{\mathrm{var}}=\displaystyle \frac{1}{{F}_{\mathrm{var}}}{\left(\displaystyle \frac{1}{2\,N}\right)}^{1/2}\displaystyle \frac{{S}^{2}}{{\bar{F}}^{2}}.\end{eqnarray} \tag{ 4 }$$

The variability characteristics of our light curves are given in Table 4. The variability amplitudes of the emission lines for both PG 0923+201 and PG 1001+291 can be ranked in the following order: Hγ > Hβ, Hγ > Hβ > Fe ii. Such variability ranks are commonly seen in previous observations (e.g., Bentz et al. 2010; Barth et al. 2015; Hu et al. 2020a, 2020b). The mean fluxes of combined continuum and broad emission lines are also listed in Table 4.

Table 4. Light-curve Statistics

Object	Light Curve	Mean Flux	Fvar(%)	${R}_{\max }$	τd	σd
PG 0923+201	Con	28.06 ± 5.12	17.81 ± 0.59	3.20	965 ± 303	3.47 ± 0.56
	Hβ	21.78 ± 2.14	9.74 ± 0.75	1.40	315 ± 169	1.88 ± 0.47
	Hγ	5.84 ± 0.73	12.06 ± 0.99	1.61	277 ± 167	0.62 ± 0.16
PG 1001+291	Con	14.90 ± 0.90	5.97 ± 0.27	1.30	279 ± 160	0.72 ± 0.19
	Hβ	8.85 ± 0.37	4.06 ± 0.23	1.28	135 ± 92	0.35 ± 0.10
	Hγ	2.46 ± 0.19	6.72 ± 0.48	1.51	158 ± 107	0.15 ± 0.04
	Fe ii	11.58 ± 0.44	3.17 ± 0.24	1.24	57 ± 38	0.37 ± 0.08

Note. The continuum fluxes are in units of 10⁻¹⁶ erg s⁻¹ cm⁻² Å⁻¹ and the fluxes of all emission lines are in units of 10⁻¹⁴ erg s⁻¹ cm⁻². Fluxes are corrected with the Galactic extinction. τ_d is the rest-frame damping timescale in units of day and σ_d is the variation amplitude, which has the same unit as the continuum/emission-line fluxes (see Section 4.2).

Download table as:  ASCIITypeset image

We calculate time lags between the combined continuum and broad emission-line flux variations through three methods: the interpolated cross-correlation function (ICCF; Gaskell & Sparke 1986; Gaskell & Peterson 1987), JAVELIN (Zu et al.2011), and MICA ²¹ (Li et al. 2016). For the ICCF method, the time lag is estimated by τ_cent, defined as the centroid of the ICCF above 80% of the peak value (${r}_{\max };$ Peterson et al. 2004). The uncertainties of the time lags are obtained by the 15.87% and 84.13% quantiles of the cross-correlation centroid distributions (CCCDs), generated by the "flux randomization/random subset sampling (FR/RSS)" method (Peterson et al. 1998).

Both JAVELIN and MICA use the damped random walk (DRW) model (e.g., Kelly et al. 2009) to describe the continuum variability and a specific transfer function to fit the light curves of emission lines. JAVELIN adopts a top-hat function to approximate the realistic transfer function, while MICA adopts a family of displaced Gaussians. For simplicity, we use only one Gaussian in MICA. We assign time lags as the centers of the top hat/Gaussian. In MICA, we additionally switch on the functionality of including a parameter for any unknown systematic errors, which is added to the data errors in quadrature. JAVELIN employs the affine invariant sampling algorithm (Goodman & Weare 2010; implemented by the package emcee ²² ) and MICA employs the diffusive nested sampling algorithm (Brewer et al. 2011; implemented by the package cdnest ²³ ) to perform the Markov Chain Monte Carlo (MCMC) technique. This generates posterior samples of the model parameters. The time lags and their associated errors are estimated from the median, and 15.87% and 84.13% quantiles of the corresponding posterior distributions. Below, we denote time lags obtained by JAVELIN and MICA as τ_JAV and τ_MICA, respectively.

In Figures 5 and 6, the left panels show the reconstructed light curves by JAVELIN (in orange) and MICA (in blue). The right panels show the autocorrelation functions (ACFs) of the combined continuum light curves, ICCFs, as well as the time lag distributions in the observed frame.

The time lags measured by the three methods are listed in Table 5. We can find general agreements to within uncertainties among the results of the three methods. In particular for PG 0923+201, the time lags of Hβ and Hγ measured from ICCF, JAVELIN, and MICA are well consistent with each other. For PG 1001+291, the Hβ, Hγ, and Fe ii lags obtained by JAVELIN and MICA are in good agreement.

Table 5. Rest-frame Time Lag and Transfer Function Width Measurements

Object	Line	${r}_{\max }$	τcent	τMICA	τJAV	WMICA	WJAV
			(day)			(day)
PG 0923+201	Hβ	0.94	${108.2}_{-12.3}^{+6.6}$	${105.4}_{-4.9}^{+4.9}$	${104.3}_{-1.4}^{+2.5}$	${40.2}_{-27.9}^{+41.6}\,({45.8}_{-33.6}^{+35.9})$	${4.9}_{-4.0}^{+15.4}\,({11.6}_{-10.8}^{+8.6})$
	Hγ	0.90	${121.0}_{-10.2}^{+8.5}$	${111.2}_{-6.5}^{+7.5}$	${113.9}_{-5.3}^{+5.3}$	${21.8}_{-11.7}^{+19.9}\,({25.5}_{-15.4}^{+16.1})$	${52.5}_{-21.1}^{+15.7}\,({50.2}_{-18.8}^{+18.0})$
PG 1001+291	Hβ	0.89	${37.3}_{-6.0}^{+6.9}$	${47.5}_{-5.1}^{+4.7}$	${48.2}_{-3.1}^{+3.3}$	${68.0}_{-9.8}^{+10.4}\,({68.3}_{-10.1}^{+10.1})$	${89.1}_{-8.0}^{+8.2}\,({89.2}_{-8.1}^{+8.1})$
	Hγ	0.71	${28.0}_{-17.6}^{+17.9}$	${55.1}_{-10.3}^{+10.3}$	${57.4}_{-10.2}^{+7.3}$	${166.9}_{-27.0}^{+29.2}\,({168.4}_{-28.5}^{+27.7})$	${165.8}_{-34.4}^{+26.2}\,({163.0}_{-31.6}^{+29.0})$
	Fe ii	0.78	${57.2}_{-11.8}^{+10.5}$	${74.5}_{-9.3}^{+8.4}$	${81.8}_{-5.5}^{+5.0}$	${115.7}_{-27.7}^{+26.6}\,({114.4}_{-26.4}^{+27.9})$	${165.8}_{-21.1}^{+19.3}\,({164.9}_{-20.2}^{+20.2})$

Note. "W_MICA" is the FWHM of the Gaussian transfer function obtained in the MICA fits, and "W_JAV" is the width of the top-hat transfer function obtained in the JAVELIN fits. The widths and uncertainties are estimated from the median and 68.3% confidence levels of the corresponding posterior distributions, respectively. The values in brackets refer to the means with the 68.3% confidence levels.

Download table as:  ASCIITypeset image

It is worth mentioning that PG 1001+291 was previously monitored by Du et al. (2018b) between 2015 and 2017, who reported an Hβ time lag of ${32.2}_{-4.2}^{+43.5}$ days and a mean 5100 Å luminosity of (3.31 ± 0.08) × 10⁴⁵ erg s⁻¹. The large upper error of the reported Hβ time lag was caused by the relatively short temporal baseline and low variation amplitudes of the Hβ light curve. By comparison, our campaign obtains a Hβ time lag of ${37.3}_{-6.0}^{+6.9}$ days and a slightly higher mean 5100 Å luminosity of (3.90 ± 0.24) × 10⁴⁵ erg s⁻¹ (due to variability). The time lag uncertainty of our measurement is more symmetric and significantly reduced.

Figure 7 shows the relationship between time lag ratios of Fe ii to Hβ (τ_{Fe II} /τ_Hβ ) and flux ratios of Fe ii to Hβ for PG 1001+291 together with the samples of Hu et al. (2015) and 3C 273 from Zhang et al. (2019). The results for PG 1001+291 are generally consistent with the trend that when ${{ \mathcal R }}_{\mathrm{Fe}}\gt 1$, τ_{Fe II} is approximately equal to τ_Hβ , whereas when ${{ \mathcal R }}_{\mathrm{Fe}}\lt 1$, τ_{Fe II} is larger than τ_Hβ .

Zoom In Zoom Out Reset image size

Through the Monte Carlo simulation tests described below, we demonstrate that our measured time lags are reliable and not caused by seasonal gaps. We adopt the τ_cent values to calculate the black hole masses in the following analysis.

Because our light curves have seasonal gaps, we employ Monte Carlo simulations to test whether the correlations between light curves are caused by seasonal gaps. We generate uncorrelated mock light curves based on the observed light curves of the continuum and broad emission line. The mock light curves follow the DRW process, which has a covariance between times t_i and t_j (e.g., Kelly et al. 2009):

$$\begin{eqnarray}&&S({t}_{i},{t}_{j})={\sigma }_{{\rm{d}}}^{2}\,\exp \left(-\displaystyle \frac{| {t}_{i}-{t}_{j}| }{{\tau }_{{\rm{d}}}}\right),\end{eqnarray} \tag{ 5 }$$

where τ_d is a damping timescale and σ_d is the variation amplitude. We first use the observed light curves to determine the parameters τ_d and σ_d and their uncertainties (listed in Table 4), from which we randomly draw pairs of τ_d and σ_d. We then generate two sets of mock light curves with daily samplings and use a linear interpolation onto the observed epochs to mimic real observations. We add Gaussian noises to the mock light curves by enforcing the relative errors equal to those of the observed light curves. We finally use the same ICCF method to calculate ${r}_{\max }$ of the two sets of light curves. This process is repeated 10,000 times, and we count the probability that ${r}_{\max }$ is higher than that of the observed light curves. We call this the false-alarm probability. For PG 0923+201, the false-alarm probabilities of Hβ and Hγ are 0.0001 and 0.0005, respectively. For PG 1001+291, the false-alarm probabilities of Hβ, Hγ, and Fe ii are 0.0007, 0.0260, and 0.0023, respectively. These quantities are relatively low, indicating that the correlations of our light curves are realistic and the influence of seasonal gaps is minor. As an example, we show the ${r}_{\max }$ distributions of mock light curves for Hβ in the left panels of Figure 8.

Zoom In Zoom Out Reset image size

In addition to the above tests, we also design tests to check whether the obtained time lags are reliable. We use the same method above to construct a daily sampled continuum light curve and convolve it with a Gaussian transfer function to generate a mock light curve of the emission line. The center and width of the Gaussian transfer function (listed in Table 5) are randomly assigned according to the best estimates and uncertainties obtained by MICA on the observed light curves. The mock light curves are again interpolated onto observed epochs. We repeat this process 10,000 times and perform the ICCF method to calculate distributions of τ_cent. Figure 8 shows an example of the distributions of Hβ time lags for the mock light curves, which are well consistent with the input values. The results for the other broad lines are similar, implying that our measured time lags are reliable.

In this section, we calculate velocity-resolved time lags for PG 0923+201 and PG 1001+291, which deliver information about the geometry and kinematics of the BLR.

With the spectral decompositions, we extract the broad Hβ component of the model fitted to each individual spectrum and calculate the rms spectrum. We then divide the rms spectrum into several velocity bins, each bin with the same integrated flux. The light curve of each bin is cross-correlated against the combined continuum to obtain time lags. The results are shown in Figure 9, where the top panels show the velocity-resolved delays and the bottom panels show the rms spectra of the broad Hβ. In Appendix B, we also show the light curves and ICCF analysis of each velocity bin. Below we comment on each object individually.

Zoom In Zoom Out Reset image size

PG 0923 + 201—The result shows shorter lags at higher velocities on both red and blue wings of the line profile, as would be expected for virialized motions. However, the lags in the redshifted bins are overall longer than those in the blueshifted bins, which might indicate a signature of outflow (Bentz et al. 2009; Du et al. 2018a). The virial envelopes in Figure 9 clearly show the asymmetric lag structure. We note that the center wavelengths of the narrow lines in our spectra (see Figure 3) are correct, meaning that this asymmetry is real. A similar pattern can be seen in other objects, such as NGC 3227 (Denney et al. 2009) and Mrk 79 (Lu et al. 2019). This complicated structure may imply the coexistence of virialization and outflow in the BLR of PG 0923+201.

PG 1001 + 291—The lag structure of PG 1001+291 is relatively symmetric and also shows smaller delays in both wings of the line profile, basically agreeing with the virial envelopes. However, this object has a low variability in each velocity bin of the Hβ line, resulting in broad ICCFs and large lag errors. Further spectroscopic monitoring on PG 1001+291 is warranted to strengthen the evidence on the velocity-delay structure.

We also use the Monte Carlo simulation method described in Section 4.2 to calculate the false-alarm probability (in terms of ${r}_{\max }$) for each wavelength bin in the delay maps. The resulting false-alarm probabilities range between 0.0001 and 0.0197 for PG 0923+201, and between 0.0022 and 0.0116 for PG 1001+291. Again, these quantities are relatively low, indicating that the correlations over wavelength bins are not dominated by seasonal gaps.

We use Equation (1) to calculate the black hole mass. There are two line-width measures, namely, FWHM and line dispersion (σ_line) of mean or rms spectra. To measure Hβ widths from the mean spectrum, we need to subtract the narrow-line components. We assume that the narrow Hβ has the same profile and a fixed flux ratio of 0.1 to [O iii] λ5007. Regarding uncertainties of the line width, we consider the following factors. First, the contribution from the narrow line is estimated by setting a flux ratio of 0 and 0.2 between the narrow Hβ to [O iii] λ5007. As such, we derive an upper and lower line width and assign the uncertainty as the average of the differences to the fiducial value with a flux ratio of 0.1. Second, we use the bootstrap method to estimate the uncertainties of line widths caused by data sampling (see, e.g., Peterson et al. 2004). Specifically, we randomly select N spectra (with replacement) from our N total spectra and remove duplicate spectra in creating a new mean spectrum. We apply the above procedure to this newly generated mean spectrum to calculate its FWHM and σ_line. We repeat this process 1000 times to obtain line-width distributions, from which we determine the standard deviations and assign them as the uncertainties. Finally, for PG 0923+201, we additionally include the influence of He ii by assigning the contributed uncertainty as the difference between the line widths with and without adding a He ii component in the spectral decomposition. We combine the above uncertainties in quadrature to get the final line-width uncertainties.

For the rms spectrum, we determine the line widths by spectral fitting, as shown in the bottom panels of Figure 3. The fitting components for the rms spectrum in PG 0923+201 include (1) a single power law, (2) a single Gaussian for Hβ, He ii, and Hγ; and in PG 1001+291 include (1) a single power law, (2) a Fe ii template from Boroson & Green (1992), and (3) a single Gaussian for Hβ and Hγ. The corresponding uncertainties are determined by the same bootstrap method as applied for the mean spectrum described above.

From all measured line widths we subtract in quadrature the instrumental broadening (Section 2.2) to obtain the results in Table 6. Following Du et al. (2018b), we adopt τ_cent and the Hβ FWHM from the mean spectrum and f_BLR = 1 to measure the black hole mass. The mass errors simply include the uncertainties of time lags and line widths. We obtain a black hole mass of ${118}_{-16}^{+11}\times {10}^{7}{M}_{\odot }$ for PG 0923+201 and ${3.33}_{-0.54}^{+0.62}\times {10}^{7}{M}_{\odot }$ for PG 1001+291.

Table 6. Hβ Line-width Measurements and Black Hole Mass Estimates

	Mean Spectrum		rms Spectrum
Object	FWHM	σline	FWHM	σline	M•	L5100	$\dot{{ \mathcal M }}$
	(km s−1)		(km s−1)		(107 M⊙)	(1045 erg s−1)
PG 0923+201	7469 ± 273	3156 ± 109	6853 ± 214	2914 ± 88	${118}_{-16}^{+11}$	2.05 ± 0.29	${0.21}_{-0.07}^{+0.06}$
PG 1001+291	2138 ± 9	1320 ± 4	2108 ± 120	896 ± 51	${3.33}_{-0.54}^{+0.62}$	3.90 ± 0.24	${679}_{-227}^{+259}$

Download table as:  ASCIITypeset image

We also measure the black hole mass using the broad Hγ line. Because the velocity width and shift of Hγ are fixed to those of Hβ in fitting the mean spectrum, we resort to the FWHM from the rms spectrum. For PG 0923+201 and PG 1001+291, the FWHMs of Hγ are 7572 ± 404 and 1884 ± 250 km s⁻¹, respectively. Again, using the width and τ_cent of Hγ and f_BLR = 1, we estimate the black hole mass to be ${135}_{-18}^{+17}\times {10}^{7}{M}_{\odot }$ for PG 0923+201 and ${1.94}_{-1.32}^{+1.34}\times {10}^{7}{M}_{\odot }$ for PG 1001+291, respectively. The obtained black hole masses from the Hβ and Hγ are consistent with each other.

According to the standard accretion disk model (Shakura & Sunyaev 1973), the dimensionless accretion rate (Eddington ratio) defined as $\dot{{\mathscr{M}}\,}\,={\dot{M}}_{\bullet }\,{c}^{2}/{L}_{{\rm{Edd}}}$ can be estimated from Du et al. (2016):

$$\begin{eqnarray}&&\dot{{\mathscr{M}}\,}\,=20.1{\left(\displaystyle \frac{{\ell }_{44}}{\cos i}\right)}^{3/2}{m}_{7}^{-2},\end{eqnarray} \tag{ 6 }$$

where ${\dot{M}}_{\bullet }$ is the accretion rate, L_Edd is the Eddington luminosity, ℓ₄₄ = L₅₁₀₀/10⁴⁴erg s⁻¹, m₇ = M_•/10⁷ M_⊙, and i is the inclination angle of the accretion disk. We take $\cos i=0.75$, which represents a mean disk inclination for a type 1 AGN by assuming that the inclination is randomly distributed between 0° and 60°. We estimate the flux contributions from the host galaxies based on the empirical relation proposed by Shen et al. (2011), which is written as L_5100,host/L_5100,AGN = 0.8052 − 1.5502x + 0.912x² − 0.1577x³, for x < 1.053, where L_5100,host and L_5100,AGN are the luminosity of the host and AGN at 5100 Å, respectively, $x=\mathrm{log}\left({L}_{5100,\mathrm{tot}}/{10}^{44}\mathrm{erg}\,{{\rm{s}}}^{-1}\right)$, and L_5100,tot is the total luminosity at 5100 Å. For x > 1.053, the luminosity correction is unnecessary. For PG 0923+201 and PG 1001+291, x = 1.380 and 1.623, respectively, meaning that the contribution of host galaxies can be ignored. This also demonstrates that it is reasonable to omit the host galaxy component in our spectral decompositions in Section 3.1. We calculate the mean flux of the continuum light curve obtained by the spectral fitting and use it to compute the luminosity at 5100 Å. Combining the M_• and L₅₁₀₀, we estimate $\dot{{\mathscr{M}}\,}$ to be ${0.21}_{-0.07}^{+0.06}$ for PG 0923+201 and ${679}_{-227}^{+259}$ for PG 1001+291. This implies that PG 0923+201 is a sub-Eddington accretor whereas PG 1001+291 is a super-Eddington accretor.

In the above calculations, we adopt a virial factor of f_BLR = 1. Because the two objects have significantly different properties, their virial factors might be different. Several previous studies proposed that the virial factor f_BLR for Hβ line is anticorrelated with the FWHM (e.g., Mejía-Restrepo et al. 2018; Martínez-Aldama et al. 2019; Yu et al. 2019). Using the relation of Mejía-Restrepo et al. (2018),

$$\begin{eqnarray}&&{f}_{\mathrm{BLR}}^{c}={\left(\displaystyle \frac{\mathrm{FWHM}}{4550\pm 1000}\right)}^{-1.17},\end{eqnarray} \tag{ 7 }$$

we obtain the virial factor based on the FWHM from the mean spectrum 0.56 ± 0.14 for PG 0923+201, and 2.42 ± 0.62 for PG 1001+291. With this new viral factor, the black hole masses are estimated to be ${66}_{-19}^{+18}\times {10}^{7}{M}_{\odot }$ for PG 0923+201 and ${8.05}_{-2.44}^{+2.55}\times {10}^{7}{M}_{\odot }$ for PG 1001+291. Correspondingly, the dimensionless accretion rates are changed to ${0.66}_{-0.40}^{+0.38}$ and ${116}_{-71}^{+74}$, respectively. This still indicates that PG 0923+201 is a sub-Eddington accretor and PG 1001+291 is a super-Eddington accretor, retaining the conclusion with f_BLR = 1.

In the left panel of Figure 10, we plot the Hβ lag and L₅₁₀₀ of PG 0923+201 and PG 1001+291 along with the previous compiled sample of Bentz et al. (2013) and samples from the SEAMBH campaigns (Du et al. 2015, 2016, 2018b). There are seven RM AGNs with luminosities above 10⁴⁵ erg s⁻¹ to date. For comparison, we also superimpose the empirical R_Hβ –L₅₁₀₀ relation from Du et al. (2018b). The dotted line represents the relation for the RM sample with $\dot{{\mathscr{M}}\,}\,\lt 3$, which is consistent with Bentz et al.'s (2013) compilation. The dashed line represents the relation for the RM sample with $\dot{{\mathscr{M}}\,}\,\geqslant 3$. The location of PG 1001+291 is almost unchanged compared to the previous measurement of Du et al. (2018b). The Hβ time delay of PG 0923+201 is consistent with the empirical R_Hβ –L₅₁₀₀ relations of both Bentz et al. (2013) and Du et al. (2018b), in consideration of their associated scatters. However, PG 1001+291 lies 0.78 dex below the empirical R_Hβ –L₅₁₀₀ relation of Bentz et al. (2013). The reported scatter σ of the Bentz et al. (2013) relation is about 0.19, and the deviation of PG 1001+291 exceeds a 4σ significance. This deviation is believed to be caused by the physical dependence of the R_Hβ –L₅₁₀₀ relation on the dimensionless accretion rates. According to the results of Du et al. (2015, 2016, 2018b), there is a strong correlation between time lag shortening in SEAMBHs and accretion rates. A possible explanation for time lag shortening was proposed by Wang et al. (2014b). Specifically, slim accretion disks with super-Eddington accretion rates produce strong self-shadowing effects, resulting in a strongly anisotropic radiation field and two dynamically distinct BLR regions with different delays (Wang et al. 2014b). The BLR is thereby divided into a shadowed region and an unshadowed region (Du et al. 2018b). Because the shadowed region receives fewer ionizing photons, its size shrinks, leading to a shortened time delay. Another possible explanation was that the time lag shortening is caused by the changes of the UV/optical spectral energy distribution and the relative amount of ionizing radiation, which may be related to the black hole spin and accretion rate (e.g., Wang et al. 2014a; Czerny et al. 2019; Fonseca Alvarez et al. 2020).

Zoom In Zoom Out Reset image size

According to the R_Hβ –L₅₁₀₀ relation of Du et al. (2018b),

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{\rm{l}}{\rm{o}}{\rm{g}}\displaystyle \frac{{R}_{{\rm{H}}\beta }}{{\rm{l}}{\rm{t}}{\rm{d}}}\\ =\,\left\{\begin{array}{ll}(1.53\pm 0.03)+(0.51\pm 0.03){\rm{l}}{\rm{o}}{\rm{g}}{\ell }_{44}, & (\dot{{\mathscr{M}}\,}\,\lt 3)\\ (1.26\pm 0.04)+(0.45\pm 0.05){\rm{l}}{\rm{o}}{\rm{g}}{\ell }_{44}, & (\dot{{\mathscr{M}}\,}\,\geqslant 3)\end{array}\right.,\end{array}\end{array}\end{eqnarray} \tag{ 8 }$$

we can find that PG 0923+201 and PG 1001+291 are located within 2σ of the relation with $\dot{{\mathscr{M}}\,}\,\lt 3$ and 3σ of the relation with $\dot{{\mathscr{M}}\,}\,\geqslant 3$, respectively. This confirms that the R_Hβ –L₅₁₀₀ relation at the high-luminosity end still strongly depends on the dimensionless accretion rates. Following Du et al. (2015), we define ${\rm{\Delta }}{R}_{{\rm{H}}\beta }=\mathrm{log}({R}_{{\rm{H}}\beta }/{R}_{{\rm{H}}\beta ,{\rm{R}}-{\rm{L}}})$ to measure the deviation from the R_Hβ –L₅₁₀₀ relation for $\dot{{\mathscr{M}}\,}\,\lt 3$ and plot ΔR_Hβ as a function of the dimensionless accretion rate $\dot{{\mathscr{M}}\,}$in the right panel of Figure 10. PG 1001+291 clearly deviates from the canonical relation by about 0.8 dex. Figure 10 demonstrates that Hβ time lags are more severely shortened for higher accretion rates. Therefore, it is important to take accretion rates into account when estimating the black hole mass of luminous AGNs at super-Eddington accretion rates.

In Figure 11, we plot the distribution of RM samples in the Eigenvector 1 (EV1) plane (also known as the main sequence) by including the compilation of Du & Wang (2019) and two objects obtained in this work. Here, the EV1 plane refers to the ${{ \mathcal R }}_{\mathrm{Fe}}$ versus FWHM_Hβ (FWHM of the broad Hβ) plane. Sulentic et al. (2000) divides AGNs into populations A and B according to the value of FWHM_Hβ , where AGNs with FWHM_Hβ ≤ 4000 km s⁻¹ are classified as Population A, otherwise as Population B. Population A includes NLS1s and high accretors (Marziani & Sulentic 2014), while Population B has larger black hole masses and lower accretion rates (Sulentic et al. 2011). It is generally believed that the accretion rate and Eddington ratio increase with ${{ \mathcal R }}_{\mathrm{Fe}}$, and the dispersion of FWHM_Hβ for a given ${{ \mathcal R }}_{\mathrm{Fe}}$ characterizes the orientation effect (Marziani et al. 2001; Shen & Ho 2014). PG 1001+291 is located in the region of high ${{ \mathcal R }}_{\mathrm{Fe}}$ and narrow FWHM_Hβ , consistent with our results that PG 1001+291 is accreting at a super-Eddington rate. We note that the ${{ \mathcal R }}_{\mathrm{Fe}}$ of PG 1001+291 is larger than the previous measurement of Du & Wang (2019) due to variability.

Zoom In Zoom Out Reset image size

We simply follow Sulentic et al. (2011) to define outliers in the EV1 plane, namely, those AGNs with FWHM_Hβ > 4000 km s⁻¹ and ${{ \mathcal R }}_{\mathrm{Fe}}\,\gt 0.5$. PG 0923+201 is the most significant outlier in the present RM sample. The reasons for such an "outlier" location of PG 0923+201 in the EV1 plane are not yet clear. However, we note that this definition of "outlier" is only phenomenological. Meanwhile, we cannot exclude the possibility that the significant deviation might be caused by real differences between PG 0923+201 and the RM sample of Du & Wang (2019). The majority of the Du & Wang (2019) sample have low luminosities and black hole masses, and therefore, the corresponding FWHMs are generally small. Nevertheless, from the large SDSS quasar sample, there also exists a population of AGNs that have relatively large ${{ \mathcal R }}_{\mathrm{Fe}}$ and broad Hβ line widths like PG 0923+201 (see Shen & Ho 2014). Detailed multiwavelength investigations of these quasars might help reveal additional BLR and SMBH accretion physics.