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

Sha-Sha Li  Key Laboratory for Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, 19B Yuquan Road, Beijing 100049,
People’s Republic of China  University of Chinese Academy of Sciences, 19A
Yuquan Road, Beijing 100049, People’s Republic of China  Sen Yang  Key
Laboratory for Particle Astrophysics, Institute of High Energy Physics,
Chinese Academy of Sciences, 19B Yuquan Road, Beijing 100049, People’s
Republic of China  University of Chinese Academy of Sciences, 19A Yuquan Road,
Beijing 100049, People’s Republic of China  Zi-Xu Yang  Key Laboratory for
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Beijing 100049, People’s Republic of China
University of Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049,
People’s Republic of China  Yong-Jie Chen  Key Laboratory for Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, People’s Republic of China  University of
Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s
Republic of China  Yu-Yang Songsheng Key Laboratory for Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan
Road, Beijing 100049, People’s Republic of China  University of Chinese
Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s Republic of
China  He-Zhen Liu  School of Astronomy and Space Science, Nanjing University,
Nanjing, Jiangsu 210093, People’s Republic of China  Key Laboratory of Modern
Astronomy and Astrophysics(Nanjing University), Ministry of Education,
Nanjing, Jiangsu 210093, People’s Republic of China  Pu Du Key Laboratory for
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Beijing 100049, People’s Republic of China  Bin Luo
School of Astronomy and Space Science, Nanjing University, Nanjing, Jiangsu
210093, People’s Republic of China  Key Laboratory of Modern Astronomy and
Astrophysics(Nanjing University), Ministry of Education, Nanjing, Jiangsu
210093, People’s Republic of China  Zhe Yu  Key Laboratory for Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, People’s Republic of China  University of
Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s
Republic of China  Chen Hu  Key Laboratory for Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan
Road, Beijing 100049, People’s Republic of China  Bo-Wei Jiang Key Laboratory
for Particle Astrophysics, Institute of High Energy Physics, Chinese Academy
of Sciences, 19B Yuquan Road, Beijing 100049, People’s Republic of China
University of Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049,
People’s Republic of China  Dong-Wei Bao  Key Laboratory for Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, People’s Republic of China  University of
Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s
Republic of China  Wei-Jian Guo  Key Laboratory for Particle Astrophysics,
Institute of High Energy Physics, Chinese Academy of Sciences, 19B Yuquan
Road, Beijing 100049, People’s Republic of China  University of Chinese
Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s Republic of
China  Zhi-Xiang Zhang  Department of Astronomy, Xiamen University, Xiamen,
Fujian 361005, People’s Republic of China  Yan-Rong Li Key Laboratory for
Particle Astrophysics, Institute of High Energy Physics, Chinese Academy of
Sciences, 19B Yuquan Road, Beijing 100049, People’s Republic of China  Ming
Xiao  Key Laboratory for Particle Astrophysics, Institute of High Energy
Physics, Chinese Academy of Sciences, 19B Yuquan Road, Beijing 100049,
People’s Republic of China  Kai-Xing Lu Yunnan Observatories, Chinese Academy
of Sciences, Kunming 650011, People’s Republic of China  Luis C. Ho Kavli
Institute for Astronomy and Astrophysics, Peking University, Beijing 100871,
People’s Republic of China  Department of Astronomy, School of Physics, Peking
University, Beijing 100871, People’s Republic of China  Jing-Min Bai  Yunnan
Observatories, Chinese Academy of Sciences, Kunming 650011, People’s Republic
of China  Wei-Hao Bian Physics Department, Nanjing Normal University, Nanjing
210097, People’s Republic of China  Jesús Aceituno  Centro Astronomico Hispano
Alemán, Sierra de los filabres sn, E-04550 GéTablesrgal, Almería, Spain
Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la astronomía sn,
E-18008 Granada, Spain  Takeo Minezaki Institute of Astronomy, School of
Science, University of Tokyo, Mitaka, Tokyo 181-0015, Japan  Keith Horne SUPA
Physics and Astronomy, University of St. Andrews, North Haugh, KY16 9SS, UK
Mitsuru Kokubo Department of Astrophysical Sciences, Princeton University,
Princeton, New Jersey 08544,USA  Jian-Min Wang Key Laboratory for Particle
Astrophysics, Institute of High Energy Physics, Chinese Academy of Sciences,
19B Yuquan Road, Beijing 100049, People’s Republic of China  University of
Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, People’s
Republic of China  National Astronomical Observatories of China, Chinese
Academy of Sciences, 20A Datun Road, Beijing 100020, People’s Republic of
China Jian-Min Wang<EMAIL_ADDRESS>
###### Abstract

We report the results of a multi-year spectroscopic and photometric monitoring
campaign of two luminous quasars, PG 0923+201 and PG 1001+291, both located at
the high-luminosity end of the broad-line region (BLR) size-luminosity
relation with optical luminosities above $10^{45}~{}{\rm erg~{}s^{-1}}$. PG
0923+201 is for the first time monitored, and PG 1001+291 was previously
monitored but our campaign has a much longer temporal baseline. We detect time
lags of variations of the broad H$\beta$, H$\gamma$, Fe ii lines with respect
to those of the 5100 Å continuum. The velocity-resolved delay map of H$\beta$
in PG 0923+201 indicates a complicated structure with a mix of Keplerian disk-
like motion and outflow, and the map of H$\beta$ in PG 1001+291 shows a
signature of Keplerian disk-like motion. Assuming a virial factor of $f_{\rm
BLR}=1$ and FWHM line widths, we measure the black hole mass to be
$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. Their respective
accretion rates are estimated to be $0.21_{-0.07}^{+0.06}\times L_{\rm
Edd}\,c^{-2}$ and $679_{-227}^{+259}\times L_{\rm Edd}\,c^{-2}$, indicating
that PG 0923+201 is a sub-Eddington accretor and PG 1001+291 is a super-
Eddington accretor. While the H$\beta$ time lag of PG 0923+201 agrees with the
size-luminosity relation, the time lag of PG 1001+291 shows a significant
deviation, confirming that in high-luminosity AGN the BLR size depends on both
luminosity and Eddington ratio. Black hole mass estimates from single AGN
spectra will be over-estimated at high luminosities and redshifts if this
effect is not taken into account.

Quasars(1319); Supermassive black holes (1663); Reverberation mapping (2019)

##  1 Introduction

Over the past four decades, reverberation mapping (RM) technique (Bahcall et
al., 1972; Blandford & McKee, 1982; Peterson, 1993, 2014) has been established
as a standard tool for measuring the mass of the central supermassive black
hole (SMBH) and for studying geometry and dynamics of the broad-line region
(BLR) in active galactic nuclei (AGNs). The characteristic time delay
$\tau_{\rm BLR}$ between broad emission line and continuum variations
corresponds to the light traveling time from the central continuum source to
the BLR, and therefore implies the BLR size by multiplying the time delay with
the light speed, namely, $R_{\rm BLR}=c\,\tau_{\rm BLR}$. Assuming the BLR is
virialized, the black hole mass can be determined through $R_{\rm BLR}$ and
the broad emission-line width $\Delta V$,

$M_{\bullet}=f_{\rm BLR}\frac{R_{\rm BLR}\Delta V^{2}}{G},$ (1)

where $G$ is the gravitational constant and $f_{\rm BLR}$ is the so-called
virial factor depending on geometry and kinematics of the BLR (e.g., Ho & Kim
2014).

Up to now, there are about a hundred AGN RM observations and black hole mass
measurements (e.g., Peterson et al. 1998, 2002, 2004; Kaspi et al. 2000; Bentz
et al. 2008, 2009; Denney et al. 2009, 2010; Barth et al. 2011, 2013, 2015;
Grier et al. 2012, 2017; Du et al. 2014, 2015, 2016b, 2018a, 2018b; Shen et
al. 2016; Fausnaugh et al. 2017; De Rosa et al. 2018; Huang et al. 2019; Lu et
al. 2019; Zhang et al. 2019; Hu et al. 2020a, b, 2021). Those RM observations
have established the widely used $R_{\rm H\beta}-L_{5100}$ relation between
the size of H$\beta$ BLR ($R_{\rm H\beta}$) and the monochromatic luminosity
at 5100 Å ($L_{5100}$) (e.g., Kaspi et al. 2000; Peterson & Wandel 1999, 2000;
Bentz et al. 2013), which allows us to infer BLR sizes with single-epoch
spectra and therefore to economically estimate black hole masses for large AGN
samples. Needless to say, to reliably estimate black hole masses in single-
epoch observations, the RM sample used for building up the $R_{\rm
H\beta}-L_{5100}$ relation needs to cover a variety of AGN populations. Recent
RM campaigns indeed show that H$\beta$ time lags of super-Eddington accreting
AGNs deviate from the $R_{\rm H\beta}-L_{5100}$ relation (Du et al., 2015,
2016b, 2018b; Fonseca Alvarez et al., 2020).

On the other hand, due to the intensive time demands of RM observations, most
of the previous RM campaigns were concentrated on low-luminosity and low-
redshift sources. Time delays of those sources are relatively short (days or
weeks) and the required monitoring periods are only several months. By
contrast, RM of high-luminosity and high-redshift AGNs suffers the following
restrictions. First, optical variability of high-luminosity AGNs is generally
weak (e.g., Hook et al. 1994; Giveon et al. 1999; Vanden Berk et al. 2004;
Kelly et al. 2009), imposing difficulties on RM monitoring. Second, time
delays for higher-redshift and higher-luminosity AGN are generally longer and
monitoring periods are accordingly needed to be extended to years. Third, for
high-redshift AGNs ($z\gtrsim 0.7$), H$\beta$ emission lines are redshifted
out of optical bands and infrared RM is required. Unfortunately, there has not
been such infrared RM campaigns for H$\beta$ yet. Hitherto, there are only a
few RM sources at the high-luminosity ($L_{5100}>10^{45}~{}{\rm
erg~{}s^{-1}}$) end of the $R_{\rm H\beta}-L_{5100}$ relation. In this regard,
it is quite necessary to expand the high-luminosity RM sample.

To this end, we undertook a multi-year RM campaign on two luminous quasars, PG
0923+201 and PG 1001+291, both with an optical luminosity above
$10^{45}~{}{\rm erg~{}s^{-1}}$. The obtained data allow us to measure the time
delays of the broad H$\beta$, H$\gamma$, and Fe II lines as well as to
constrain the BLR kinematics of the broad H$\beta$ line.

The paper is organized as follows. Observations and data reduction are
described in Section 2. In Section 3, we explain in detail measurements of
light curves of continuum and broad emission lines as well as
intercalibrations of continuum light curves from our campaign and other time-
domain survey archives. In Section 4, we present the time delay analysis and
black hole mass measurements. In Section 5, we discuss the implication of our
results on the $R_{\rm H\beta}-L_{5100}$ relation and the location of PG
0923+201 in the main sequence of RM samples. Here, the main sequence refers to
the relation between full widths at half maximum (FWHMs) of broad H$\beta$
lines and strengths of Fe ii (${\cal R}_{\rm Fe}$), where ${\cal R}_{\rm Fe}$
is the flux ratio of Fe ii emission lines between 4434 Å and 4684 Å to broad
H$\beta$ lines. The conclusions are summarized in Section 6. Throughout the
paper, we use a $\Lambda$CDM cosmology with $H_{0}=67\,\rm km\,s^{-1}\,{\rm
Mpc}^{-1}$, $\Omega_{\rm M}=0.32$, and $\Omega_{\Lambda}=0.68$ (Planck
Collaboration et al., 2020).

Figure 1: Images of PG 0923+201 (left) and PG 1001+291 (right) and their
comparison stars. Stars 1-4 in the image of PG 0923+201 and stars 1-3 in the
image of PG 1001+291 are selected for differential photometry. The blue dotted
lines illustrate the directions of slits.

Figure 2: Photometric light curves of the comparison stars for PG 0923+201
(top) and PG 1001+291 (bottom), in units of instrumental magnitudes. Blue
dashed lines represent the mean magnitudes and orange dashed lines represent
the standard deviations.

Figure 3: Spectral decompositions of the mean spectrum (top) and the rms
spectrum (bottom) for PG 0923+201 in left panels and PG 1001+291 in right
panels. In each case, two vertical adjacent panels show the details of
spectral decompositions and the residuals in fitting windows. The vertical
dotted lines indicate the regions omitted from the fits. Table 1: Properties
of spectral and photometric data

Object | Source | Monitoring Period |  | Epochs | $\Delta T$ (day)
---|---|---|---|---|---
JD - 2457000 | Date |
PG 0923+201 | LJS | 1054 $-$ 1975 | 2017 Oct $-$ 2020 May |  | 89 | 5.9
| LJP | 1054 $-$ 2013 | 2017 Oct $-$ 2020 Jun |  | 85 | 6.0
| ASAS-SN | 298 $-$ 2021 | 2015 Oct $-$ 2020 Jun |  | 519 | 1.7
| ZTF | 1202 $-$ 2000 | 2018 Mar $-$ 2020 May |  | 374 | 0.1
PG 1001+291 | LJS | 334 $-$ 1632 | 2015 Nov $-$ 2019 May |  | 164 | 4.9
| LJP | 333 $-$ 1632 | 2015 Nov $-$ 2019 Jun |  | 160 | 4.1
| ZTF | 1202 $-$ 1649 | 2018 Mar $-$ 2019 Jun |  | 143 | 1.0
| CAHAS | 894 $-$ 1492 | 2017 May $-$ 2019 Jan |  | 16 | 11.9
| CAHAP | 909 $-$ 1492 | 2017 Jun $-$ 2019 Jan |  | 14 | 11.1

Note. — LJS and CAHAS refer to spectroscopy from Lijiang and CAHA. LJP, CAHAP,
ASAS-SN, and ZTF refer to photometry from Lijiang, CAHA, ASAS-SN archive, and
ZTF archive, respectively. $\Delta T$ refers to the median sampling interval.

##  2 Observations and Data Reduction

###  2.1 Targets

We spectroscopically monitored two quasars, PG 0923+201 and PG 1001+291, both
with an optical luminosity larger than $10^{45}~{}{\rm erg~{}s^{-1}}$. Their
basic properties are summarized below.

PG 0923+201 has a redshift of $z=0.1921$, $m_{\rm g}=15.6~{}\rm mag$, and
$V$-band extinction $A_{\rm V}=0.116\,\rm mag$ (Schlafly & Finkbeiner, 2011).
Its spectrum shows prominent Fe ii emission lines, and weak [O iii]
$\lambda$5007 line (Boroson & Green, 1992). The ratios of the equivalent width
of [O iii] $\lambda$5007 and Fe ii between 4434 Å and 4684 Å to that of
H$\beta$ 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$\beta$ line is as large as $\sim 7600~{}\rm
km\,s^{-1}$, much broader than NLS1s (usually narrower than $2000~{}\rm
km\,s^{-1}$).

PG 1001+291 has a redshift of $z=0.3272$, $m_{\rm g}=15.9~{}\rm mag$, and
$A_{\rm V}=0.06~{}\rm mag$ (Schlafly & Finkbeiner, 2011). Similar to PG
0923+201, PG 1001+291 also shows prominent Fe ii emission lines with ${\cal
R}_{\rm Fe}$ $\sim$ 1.17 and weak [O iii] line (Du & Wang, 2019). However, its
H$\beta$ line has a much narrower FWHM, only about $2100~{}\rm km\,s^{-1}$. PG
1001+291 was monitored between 2015 and 2017 by the Super-Eddington Accreting
Massive Black Holes (SEAMBH) campaign111The selection criteria for SEAMBHs
were: strong optical Fe ii, relatively narrow H$\beta$, and weak [O iii] lines
(see Du et al. 2014, 2018b). (Du et al. 2018b, who used the alternative SDSS
name for PG 1001+291, namely, SDSS J100402.61+285535.3). The H$\beta$ lag
($32.2_{-4.2}^{+43.5}$ days) was relatively uncertain due to the short
monitoring period. Nevertheless, the H$\beta$ lag of PG 1001+291 was found to
be 0.80 dex below the $R_{\rm H\beta}-L_{5100}$ relation (Du et al., 2018b).

###  2.2 Spectroscopy

In our observations, spectroscopy was mainly undertaken with the Lijiang 2.4 m
telescope at the Yunnan Observatories of Chinese Academy of Sciences. It is
equipped with the Yunnan Faint Object Spectrograph and Camera (YFOSC) that 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 observation setup below. (1) For PG 0923+201,
the spectra were obtained by Grism 14 with a dispersion of 1.78 Å pixel-1 and
a wavelength coverage of $3800-7400$ Å. We used a long slit with a projected
width of $2.^{\prime\prime}5$. The instrumental broadening was about 500 km
s-1 (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 the contamination. Grim 3 has a
spectral coverage of $3800-9000$ Å and a dispersion of 2.9 Å pixel-1. We
adopted a $5.^{\prime\prime}0$ wide slit and the instrumental broadening is
about 1200 km s-1 (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 $8.^{\prime\prime}5$ and a background region
of $7.^{\prime\prime}4-14^{\prime\prime}$ 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 October 2017 and May 2020. For PG 1001+291, we obtained a total of 164
epochs between November 2015 and May 2019. 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 $\sim$ 57 for PG 0923+201 and $\sim$ 48 for PG 1001+291.

Since May 2017, 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, b). 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-1 and a
spectral coverage of $4000-8500$ Å. We used a long slit with a projected width
of $3.^{\prime\prime}0$ and the resulting instrumental broadening is about
1000 km s-1 (Hu et al., 2020b). We selected the same comparison star as in our
Lijiang observations. A uniform aperture of $10.^{\prime\prime}6$ and a
background region of $13.^{\prime\prime}8-26.^{\prime\prime}5$ were used for
spectrum extraction (see Hu et al. 2020b for data reduction details). In
total, we obtained 16 epochs between May 2017 and January 2019 at the Calar
Alto Observatory. The typical S/N per pixel at 5100 Å (rest frame) is $\sim$
68 of each night.

The Fe ii blends are relatively strong beneath the [O iii] line and the [O
iii] $\lambda$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. We note that the absolute
spectrum of the comparison star is not important, because RM analysis only
depend on relative flux variations. 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.

###  2.3 Photometry

Photometry was obtained directly through the image mode of YFOSC. We used the
Johnson $V$ filter and 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 ${\rm 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 $4.^{\prime\prime}5$ and $4.^{\prime\prime}2$ provides the best photometry
for the two objects, respectively. The inner and outer radius of the
background was uniformly set to be $8.^{\prime\prime}5-17^{\prime\prime}$.
Figure 2 shows the photometric light curves of the comparison stars over our
monitoring periods, both of which have a standard deviation $\lesssim$ 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 photometric
light curves of the two targets.

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 June 2017 and January 2019 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-SN222http://www.astronomy.ohio-state.edu/asassn/index.shtml is a
long-term project designed to survey the whole visible sky every night down to
about 17 magnitude 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 thermo-electrically cooled CCD camera. The field of view of
each camera is about 4.5 $\rm deg^{2}$, and the pixel scale is
$8^{\prime\prime}$. Observations generally used $V$ or $g$ band filters for
three dithered 90s exposures. The photometric data were processed using IRAF
apphot package with an aperture radius of $16^{\prime\prime}$ 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 October 2015 to June 2020 (see Table 1). For PG 1001+291, the ASAS-SN
data are not adopted because of the relatively poor data quality.

The ZTF333https://www.ztf.caltech.edu is designed for transients and variables
and uses the Palomar 48-inch Schmidt Telescope with a 47 square degree 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
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 March 2018 and May 2020 for PG 0923+201 and 143 epochs with a median
sampling interval of 1.0 days from March 2018 to June 2019 for PG 1001+291
(see Table 1).

Table 2: Light curves of PG 0923+201

Spectra |  | Photometry
---|---|---
JD - 2450000 | $F_{\rm 5100}$ | $F_{\rm H\beta}$ | $F_{\rm H\gamma}$ | $\rm Obs.$ |  | JD - 2450000 | mag | $\rm Obs.$
8054.37 | $26.45\pm 0.37$ | $21.46\pm 0.33$ | $5.74\pm 0.19$ | L |  | 7298.12 | $15.457\pm 0.086$ | A
8065.37 | $26.48\pm 0.38$ | $21.84\pm 0.36$ | $5.71\pm 0.22$ | L |  | 7299.12 | $15.397\pm 0.078$ | A
8069.43 | $26.81\pm 0.36$ | $20.86\pm 0.31$ | $5.44\pm 0.17$ | L |  | 7308.11 | $15.416\pm 0.054$ | A
8073.35 | $27.75\pm 0.37$ | $21.12\pm 0.32$ | $5.86\pm 0.19$ | L |  | 7310.11 | $15.464\pm 0.062$ | A
8077.37 | $27.47\pm 0.37$ | $20.92\pm 0.32$ | $5.56\pm 0.19$ | L |  | 7311.11 | $15.482\pm 0.058$ | A

Note. — The 5100 Å continuum flux densities are in units of $10^{-16}~{}\rm
erg\,s^{-1}\,cm^{-2}$ Å-1, and the fluxes of emission lines are in units of
$10^{-14}~{}\rm erg\,s^{-1}\,cm^{-2}$. In the “Obs.” column, “L” refers to
Lijiang, “A” refers to ASAS-SN, and “Z” refers to ZTF.

(This table is available in its entirety in a machine-readable form in the
online journal.)

Table 3: Light curves of PG 1001+291

Spectra |  | Photometry
---|---|---
JD - 2450000 | $F_{\rm 5100}$ | $F_{\rm H\beta}$ | $F_{\rm H\gamma}$ | $F_{\rm Fe\ II}$ | $\rm Obs.$ |  | JD - 2450000 | mag | $\rm Obs.$
7334.41 | $15.32\pm 0.11$ | $8.60\pm 0.11$ | $2.36\pm 0.10$ | $11.20\pm 0.23$ | $\rm L$ |  | 7333.42 | $15.944\pm 0.003$ | $\rm L$
7339.40 | $14.78\pm 0.12$ | $8.65\pm 0.13$ | $2.20\pm 0.11$ | $11.38\pm 0.27$ | $\rm L$ |  | 7334.40 | $15.945\pm 0.003$ | $\rm L$
7343.34 | $14.70\pm 0.13$ | $8.29\pm 0.14$ | $2.34\pm 0.12$ | $11.13\pm 0.29$ | $\rm L$ |  | 7339.38 | $15.956\pm 0.003$ | $\rm L$
7349.38 | $15.11\pm 0.10$ | $8.49\pm 0.07$ | $2.41\pm 0.06$ | $10.71\pm 0.15$ | $\rm L$ |  | 7343.33 | $15.962\pm 0.003$ | $\rm L$
7351.42 | $14.96\pm 0.11$ | $8.72\pm 0.10$ | $2.31\pm 0.08$ | $11.38\pm 0.19$ | $\rm L$ |  | 7349.37 | $15.970\pm 0.002$ | $\rm L$

Note. — The 5100 Å continuum flux densities are in units of $10^{-16}~{}\rm
erg\,s^{-1}\,cm^{-2}$ Å-1, and the fluxes of emission lines are in units of
$10^{-14}~{}\rm erg\,s^{-1}\,cm^{-2}$. In the “Obs.” column, “L” refers to
Lijiang, “C” refers to CAHA, and “Z” refers to ZTF.
(This table is available in its entirety in a machine-readable form in the
online journal.)

##  3 Measurements

###  3.1 Light Curves

We measure the 5100 Å flux density and broad H$\beta$ line via a multi-
component 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
DASpec444DASpec is available at https://github.com/PuDu-Astro/DASpec. software
to perform multi-component 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 $\sim$ 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) an Fe ii template from Boroson & Green (1992), (3) a single
Gaussian for broad emission lines H$\beta$, H$\gamma$, (4) double Gaussians
for [O iii] $\lambda\lambda$4959, 5007 and [O iii] $\lambda$4363 with the same
tied velocity shifts and widths among the three line. The narrow H$\beta$ is
not added because the flux of the narrow H$\beta$ is less than 2% of the total
H$\beta$ flux, and the narrow H$\beta$ flux is estimated by assuming that a
typical flux ratio of 0.1 between the narrow H$\beta$ and [O iii]
$\lambda$5007 (e.g., Veilleux & Osterbrock 1987; Hu et al. 2015).

Because the H$\gamma$ is blended with [O iii] $\lambda$4363, we fix the
velocity width and shift of H$\gamma$ to those of H$\beta$. The flux ratio of
[O iii] $\lambda$4363 relative to [O iii] $\lambda$5007 is fixed to 0.251
obtained by the best fit of the mean spectrum. The [O iii] lines have
asymmetric blue wings, which is 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 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
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 $\lambda 4686$ line is relatively prominent in
the root mean square (rms) spectrum, it is weak and highly blended with the Fe
ii and the blue wing of the H$\beta$ 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$\beta$ and broad
H$\gamma$, (4) a single Gaussian with the same velocity width and shift for
narrow lines [O iii] $\lambda\lambda$4959, 5007 and [O iii] $\lambda$4363\.
Similar to PG 0923+201, the width and shift of H$\gamma$ are fixed to those of
H$\beta$, and the flux ratio of [O iii] $\lambda$4363 to [O iii] $\lambda$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_{\rm 5100}$) and
the broad emission lines (H$\beta$ and H$\gamma$ of PG 0923+201; H$\beta$,
H$\gamma$, 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 is not enough to detect a reliable time delay so that we do not
present Fe ii analysis for PG 0923+201 in this work.

Figure 4: The intercalibrated continuum light curves that combine the 5100 Å
continuum ($F_{\rm 5100}$), and photometry from our Lijiang observations
(JV/SDSS $r^{\prime}$) and ASAS-SN and ZTF archives. “P” and “S” refer to
photometric data and 5100 Å continuum from CAHA observations, respectively.

###  3.2 Intercalibration

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 555PyCALI is available at
https://github.com/LiyrAstroph/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 the panel (a) of Figures 5 and 6.

For PG 1001+291, we obtained spectra from both Lijiang and CAHA observations,
therefore, an intercalibration of the spectra is also required. Since the
instrumental broadening width (in terms of FWHM) of Lijiang and CAHA
observations are roughly similar and far smaller than broad line widths, there
is no need to adjust the spectral resolutions. We simply use the
intercalibration factors obtained above to align fluxes of the spectra. In the
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.

Figure 5: Left panels show the light curves of the combined continuum (a) and
broad emission lines (b)-(c). Superposed are the reconstructions by JAVELIN
(in orange) and MICA (in blue). The units of the continuum and emission line
light curves are $10^{-16}~{}\rm erg\,s^{-1}\,cm^{-2}$ Å-1 and $10^{-14}~{}\rm
erg\,s^{-1}\,cm^{-2}$, respectively. Right panels show time lag analysis.
Panel (e) plots the ACF of the combined continuum, and the CCFs between the
broad emission lines with continuum light curves. Panels (f)-(g) show the
distributions of time delays for the broad emission lines. Figure 6: Left
panels show the light curves of the combined continuum (a) and broad emission
lines (b)-(d). The black and blue points with error bars are from Lijiang and
CAHA observations, respectively. Superposed are the reconstructions by JAVELIN
(in orange) and MICA (in blue). The units of the continuum and emission line
light curves are $10^{-16}~{}\rm erg\,s^{-1}\,cm^{-2}$ Å-1 and $10^{-14}~{}\rm
erg\,s^{-1}\,cm^{-2}$, respectively. Right panels show time lag analysis.
Panel (e) plots the ACF of the combined continuum, and the CCFs between the
broad emission lines with continuum light curves. Panels (f)-(h) show the
distributions of time delays for the broad emission lines.

###  3.3 Variability Characteristics

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

$F_{\rm var}=\frac{(S^{2}-\triangle^{2})^{1/2}}{\bar{F}},$ (2)

where $S^{2}$ is the sample variance, $\triangle^{2}$ is the mean square
error, and $\bar{F}$ is the sample mean flux

$S^{2}=\sum_{i=1}^{N}\frac{(F_{i}-\bar{F})^{2}}{(N-1)},\quad\triangle^{2}=\sum_{i=1}^{N}\frac{\Delta_{i}^{2}}{N},\quad\bar{F}=\sum_{i=1}^{N}\frac{F_{i}}{N},$
(3)

where $F_{i}$ is the flux of $i$th observation, $\Delta_{i}$ is the
uncertainty of $F_{i}$, and $N$ is the total number of epochs. The standard
deviation of $F_{\rm var}$ is estimated by (Edelson et al. 2002)

$\sigma_{\rm var}=\frac{1}{F_{\rm
var}}\left(\frac{1}{2\,N}\right)^{1/2}\frac{S^{2}}{{\bar{F}}^{2}}.$ (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$\gamma$ >H$\beta$, H$\gamma$
>H$\beta$ >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, b).
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 | $F_{\rm var}(\%)$ | $R_{\rm max}$ | $\tau_{\rm d}$ | $\sigma_{\rm d}$
---|---|---|---|---|---|---
PG 0923+201 | $\rm Con$ | $28.06\pm 5.12$ | $17.81\pm 0.59$ | 3.20 | $965\pm 303$ |  $3.47\pm$0.56
| ${\rm H\beta}$ | $21.78\pm 2.14$ | $9.74\pm 0.75$ | 1.40 | $315\pm 169$ |  $1.88\pm$0.47
| ${\rm H\gamma}$ | $5.84\pm 0.73$ | $12.06\pm 0.99$ | 1.61 | $277\pm 167$ |  $0.62\pm$0.16
PG 1001+291 | $\rm Con$ | $14.90\pm 0.90$ | $5.97\pm 0.27$ | 1.30 | $279\pm 160$ |  $0.72\pm$0.19
| ${\rm H\beta}$ | $8.85\pm 0.37$ | $4.06\pm 0.23$ | 1.28 | $135\pm 92$ |  $0.35\pm$0.10
| ${\rm H\gamma}$ | $2.46\pm 0.19$ | $6.72\pm 0.48$ | 1.51 | $158\pm 107$ |  $0.15\pm$0.04
| ${\rm Fe\ II}$ | $11.58\pm 0.44$ | $3.17\pm 0.24$ | 1.24 | $57\pm 38$ |  $0.37\pm$0.08

Note. — The continuum fluxes are in a unit of $10^{-16}~{}\rm
erg\,s^{-1}\,cm^{-2}$ Å-1 and the fluxes of all emission lines are in a unit
of $10^{-14}~{}\rm erg\,s^{-1}\,cm^{-2}$. Fluxes are corrected with the
Galactic extinction. $\tau_{\rm d}$ is the rest-frame damping timescale in a
unit of day and $\sigma_{\rm d}$ is the variation amplitude, which has the
same unit as the continuum/emission line fluxes (see Section 4.2).

##  4 Analysis and Results

###  4.1 Time Lags

We calculate time lags between the combined continuum and broad emission-lines
flux variations through three methods: the interpolated cross-correlation
function (ICCF; Gaskell & Sparke, 1986; Gaskell & Peterson, 1987), JAVELIN (Zu
et al., 2011), and MICA666MICA is available at
https://github.com/LiyrAstroph/MICA2. (Li et al., 2016). For the ICCF method,
the time lag is estimated by $\tau_{\rm cent}$, defined as the centroid of the
ICCF above 80$\%$ of the peak value ($r_{\rm max}$; Peterson et al. 2004). The
uncertainties of 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 emcee777Available at https://github.com/dfm/emcee.)
and MICA employs the diffusive nested sampling algorithm (Brewer et al. 2011;
implemented by the package cdnest888Available at
https://gitub.com/LiyrAstroph/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,
$15.87\%$ and $84.13\%$ quantiles of the corresponding posterior
distributions. Below, we denote time lags obtained by JAVELIN and MICA as
$\tau_{\rm JAV}$ and $\tau_{\rm 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 lags distributions in the observed frame.

The measured time lags 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$\beta$ and
H$\gamma$ measured from ICCF, JAVELIN and MICA are well consistent with each
other. For PG 1001+291, the H$\beta$, H$\gamma$ and Fe ii lags obtained by
JAVELIN and MICA are in good agreement.

It is worth mentioning that PG 1001+291 was previously monitored by Du et al.
(2018b) between 2015 and 2017, who reported an H$\beta$ time lag of
$32.2_{-4.2}^{+43.5}$ days and a mean 5100 Å luminosity of $(3.31\pm
0.08)\times 10^{45}~{}{\rm erg~{}s^{-1}}$. The large upper error of the
reported H$\beta$ time lag was caused by the relatively short temporal
baseline and low variation amplitudes of the H$\beta$ light curve. By
comparison, our campaign obtains a H$\beta$ time lag of $37.3_{-6.0}^{+6.9}$
days and a slightly higher mean 5100 Å luminosity of $(3.90\pm 0.24)\times
10^{45}~{}{\rm erg~{}s^{-1}}$ (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$\beta$
($\rm\tau_{Fe\ II}\,/\tau_{H\beta}$) and flux ratios of Fe ii to H$\beta$ for
PG 1001+291 together with the samples of Hu et al. (2015) and 3C 273 from
Zhang et al. (2019). The results of PG 1001+291 are generally consistent with
the trend that when ${{\cal R}_{\rm Fe}}>1$, $\rm\tau_{Fe\ II}$ is
approximately equal to $\rm\tau_{H\beta}$, whereas when ${{\cal R}_{\rm
Fe}}<1$, $\rm\tau_{Fe\ II}$ is larger than $\rm\tau_{H\beta}$.

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 $\tau_{\rm cent}$ values to calculate the black hole masses in the
following analysis.

Table 5: Rest-frame time lag and transfer function width measurements

Object | $\rm Line$ | $r_{\rm max}$ |  | $\tau_{\rm cent}$ | $\tau_{\rm MICA}$ | $\tau_{\rm JAV}$ |  | $W_{\rm MICA}$ | $W_{\rm JAV}$
---|---|---|---|---|---|---|---|---|---
| (day) |  | (day)
PG 0923+201 | ${\rm H\beta}$ | 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})$
| ${\rm H\gamma}$ | 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 | ${\rm H\beta}$ | 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})$
| ${\rm H\gamma}$ | 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})$
| ${\rm 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_{\rm MICA}$” is the FWHM of the Gaussian transfer function
obtained in the MICA fits, and “$W_{\rm 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.

Figure 7: Relation between time lag ratios of Fe ii to H$\beta$ and intensity
ratios of Fe ii to H$\beta$.

###  4.2 Validity Tests of ICCF Time Lags

Since 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)

$S(t_{i},t_{j})=\sigma_{\rm d}^{2}\,\exp\left(-\frac{|t_{i}-t_{j}|}{\tau_{\rm
d}}\right),$ (5)

where $\tau_{\rm d}$ is a damping timescale and $\sigma_{\rm d}$ is variation
amplitude. We first use the observed light curves to determine the parameters
$\tau_{\rm d}$ and $\sigma_{\rm d}$ and their uncertainties (listed in Table
4), from which we randomly draw pairs of $\tau_{\rm d}$ and $\sigma_{\rm 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 these of the observed light curves. We finally use the same ICCF
method to calculate $r_{\rm max}$ of the two sets of light curves. This
process is repeated 10000 times and we count the probability that $r_{\rm
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 the
H$\beta$ and H$\gamma$ are 0.0001 and 0.0005, respectively. For PG 1001+291,
the false-alarm probabilities of the H$\beta$, H$\gamma$ 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_{\rm max}$
distributions of mock light curves for the H$\beta$ in the left panels of
Figure 8.

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 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 10000 times and perform the ICCF
method to calculate distributions of $\tau_{\rm cent}$. Figure 8 shows an
example of the distributions of H$\beta$ time lags for 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.

Figure 8: Validity tests of H$\beta$ time lags for PG 0923+201 (top) and PG
1001+291 (bottom). Left panels show $r_{\rm max}$ distributions of
uncorrelated mock light curves from 10000 simulations. The orange dotted lines
represent $r_{\rm max}$ of the observed light curves. Right panels show
$\tau_{\rm cent}$ distributions of mock light curves. The orange dotted lines
represent the input time lags (see Section 4.2 for details).

Figure 9: Velocity-resolved lags of the H$\beta$ line in PG 0923+201 (left) and PG 1001+291 (right). The top panels show delays in the observed frame across velocity bins and the bottom panels show the rms spectra of the broad H$\beta$. The vertical dashed lines show the boundaries of velocity bins. The delay in each bin is calculated by the ICCF method. The horizontal dotted lines and orange band show the mean H$\beta$ centroid lag and the associated uncertainties measured in Section 4.1. The blue dotted curves show the virial envelopes for a Keplerian disk with an inclination of $\cos i=0.75$ and the estimated black hole mass in Table 6. Table 6: H$\beta$ line width measurements and black hole mass estimates  | Mean Spectrum |  | rms Spectrum |  |  |  |
---|---|---|---|---|---|---|---
Object |  FWHM |  $\sigma_{\rm line}$ |  |  FWHM |  $\sigma_{\rm line}$ |  | $M_{\bullet}$ | $L_{5100}$ |  $\dot{\mathscr{M}}$
| (km s-1) |  | (km s-1) |  |  $(10^{7}M_{\odot})$ |  $(10^{45}\,\rm erg\,s^{-1})$ |
PG 0923+201 | $7469\pm 273$ | $3156\pm 109$ |  | $6853\pm 214$ | $2914\pm 88$ |  | $118_{-16}^{+11}$ | $2.05\pm 0.29$ | $0.21_{-0.07}^{+0.06}$
PG 1001+291 | $2138\pm 9$ | $1320\pm 4$ |  | $2108\pm 120$ | $896\pm 51$ |  | $3.33_{-0.54}^{+0.62}$ | $3.90\pm 0.24$ | $679_{-227}^{+259}$

###  4.3 Velocity-resolved Delays

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$\beta$ 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$\beta$. In Appendix, we also show the
light curves and ICCF analysis of each velocity bin. Below we comment on each
object individually.

#### 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 these 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 virilization 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$\beta$ 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_{\rm 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.

###  4.4 Black Hole Masses and Accretion Rates

We use Equation (1) to calculate the black hole mass. There are two line width
measures, namely, FWHM and line dispersion ($\sigma_{\rm line}$) of mean or
rms spectra. To measure H$\beta$ widths from the mean spectrum, we need to
subtract the narrow line components. We assume that the narrow H$\beta$ has
the same profile and a fixed flux ratio of 0.1 to the [O iii] $\lambda$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$\beta$ to [O iii] $\lambda$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 $\sigma_{\rm 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 the 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 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$\beta$, He ii, H$\gamma$; and in PG 1001+291 include: (1) a single power
law, (2) an Fe ii template from Boroson & Green (1992), (3) a single Gaussian
for H$\beta$ and H$\gamma$. 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 $\tau_{\rm cent}$ and the H$\beta$ FWHM from the mean
spectrum and $f_{\rm 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.

We also measure the black hole mass using the broad H$\gamma$ line. Because
the velocity width and shift of H$\gamma$ are fixed to those of H$\beta$ 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$\gamma$ are $7572\pm 404$ and
$1884\pm 250$ km s-1, respectively. Again, using the width and $\tau_{\rm
cent}$ of H$\gamma$ and $f_{\rm 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$\beta$ and H$\gamma$ are consistent with
each other.

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

$\dot{\mathscr{M}}=20.1\,\left(\frac{\ell_{44}}{\cos
i}\right)^{3/2}m_{7}^{-2},$ (6)

where $\dot{M}_{\bullet}$ is the accretion rate, $L_{\rm Edd}$ is the
Eddington luminosity, $\ell_{44}=L_{5100}/10^{44}{\rm\rm erg\,s^{-1}}$,
$m_{7}=M_{\bullet}/10^{7}M_{\odot}$, 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 randomly
distributes between 0 and 60 degrees. 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_{\rm 5100,host}/L_{\rm
5100,AGN}=0.8052-1.5502x+0.912x^{2}-0.1577x^{3}$, for $x<1.053$, where $L_{\rm
5100,host}$ and $L_{\rm 5100,AGN}$ are the luminosity of host and AGN at 5100
Å respectively, $x=\log\left(L_{\rm 5100,tot}/10^{44}{\rm
erg\,s^{-1}}\right)$, and $L_{\rm 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 neglected. 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_{\bullet}$ and $L_{5100}$, 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_{\rm BLR}=1$. Since
the two objects have significantly different properties, their virial factors
might be different. Several previous studies proposed that the virial factor
$f_{\rm BLR}$ for H$\beta$ line is anticorrelated with the FWHM (e.g., Mejía-
Restrepo et al. 2018; Yu et al. 2019; Martínez-Aldama et al. 2019). Using the
relation of Mejía-Restrepo et al. (2018)

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

we obtain the virial factor based on the FWHM from the mean spectrum $0.56\pm
0.14$ for PG 0923+201, and $2.42\pm 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_{\rm BLR}=1$.

Figure 10: (Left) the $R_{\rm H\beta}-L_{5100}$ relation. The dotted and
dashed lines are the best fit to the AGNs with $\dot{\mathscr{M}}<3$ and with
$\dot{\mathscr{M}}\geq 3$ from Du et al. (2018b). The RM objects with
luminosities above $10^{45}~{}{\rm erg~{}s^{-1}}$ are are located to the right
of the vertical dashed line. (Right) the relation between $\Delta R_{\rm
H\beta}$ and $\dot{\mathscr{M}}$. The vertical and horizontal dashed lines
indicate $\dot{\mathscr{M}}=3$ and $\Delta R_{\rm H\beta}=0$, respectively
(see Section 5.1 for the definition of $\Delta R_{\rm H\beta}$).

##  5 Discussion

###  5.1 Implication for $R_{\rm H\beta}-L_{5100}$ Relation

In the left panel of Figure 10, we plot the H$\beta$ lag and $L_{5100}$ 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, 2016b,
2018b). There are seven RM AGNs with luminosities above $10^{\rm 45}~{}\rm
erg\,s^{-1}$ to date. For comparison, we also superimpose the empirical
$R_{\rm H\beta}-L_{5100}$ relation from Du et al. (2018b). The dotted line
represents the relation for RM sample with $\dot{\mathscr{M}}\\!<\\!3$, which
is consistent with Bentz et al. (2013)’s compilation. The dashed line
represents the relation for RM sample with $\dot{\mathscr{M}}\\!\geq\\!3$. The
location of PG 1001+291 is almost unchanged compared to the previous
measurement of Du et al. (2018b). The H$\beta$ time delay of PG 0923+201 is
consistent with the empirical $R_{\rm H\beta}-L_{5100}$ 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_{\rm H\beta}-L_{5100}$ relation of Bentz et al. (2013). The reported
scatter $\sigma$ of the Bentz et al. (2013) relation is about 0.19 and the
deviation of PG 1001+291 exceeds a $4\sigma$ significance. This deviation is
believed to be caused by the physical dependence of the $R_{\rm
H\beta}-L_{5100}$ relation on the dimensionless accretion rates. According to
the results of Du et al. (2015, 2016b, 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). Since 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).

According to the $R_{\rm H\beta}-L_{5100}$ relation of Du et al. (2018b),

$\log\left(\frac{R_{\rm H\beta}}{{\rm
ltd}}\right)\\!=\\!\left\\{\\!\begin{array}[]{ll}\\!(1.53\pm
0.03)\\!+\\!(0.51\pm 0.03)\log\ell_{44},&\\!(\dot{\mathscr{M}}\\!<\\!3)\\\ \\\
\\!(1.26\pm 0.04)\\!+\\!(0.45\pm
0.05)\log\ell_{44},&\\!(\dot{\mathscr{M}}\\!\geq\\!3)\end{array}\right.$ (8)

we can find that PG 0923+201 and PG 1001+291 are located within $2\sigma$ of
the relation with $\dot{\mathscr{M}}\\!<\\!3$ and $3\sigma$ of the relation
with $\dot{\mathscr{M}}\\!\geq\\!3$, respectively. This confirms that the
$R_{\rm H\beta}-L_{5100}$ relation at the high-luminosity end still strongly
depends on the dimensionless accretion rates. Following Du et al. (2015), we
define $\Delta R_{\rm H\beta}=\rm log(R_{\rm H\beta}/R_{\rm{H\beta},R-L})$ to
measure the deviation from the $R_{\rm H\beta}-L_{5100}$ relation for
$\dot{\mathscr{M}}<3$ and plot $\Delta R_{\rm H\beta}$ 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$\beta$ 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.

###  5.2 Locations in the Eigenvector 1 Plane of RM Samples

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 ${\cal R}_{\rm Fe}$ versus $\rm FWHM_{H\beta}$ (FWHM of the
broad H$\beta$) plane. Sulentic et al. (2000) divides AGNs into populations A
and B according to the value of $\rm FWHM_{H\beta}$, where AGNs with $\rm
FWHM_{H\beta}\leq 4000~{}km~{}s^{-1}$ 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 ${\cal R}_{\rm Fe}$,
and the dispersion of $\rm FWHM_{H\beta}$ for a given ${\cal R}_{\rm Fe}$
characterizes the orientation effect (Marziani et al., 2001; Shen & Ho, 2014).
PG 1001+291 is located in the region of high ${\cal R}_{\rm Fe}$ and narrow
$\rm FWHM_{H\beta}$, consistent with our results that PG 1001+291 is accreting
at a super-Eddington rate. We note that the ${\cal R}_{\rm Fe}$ of PG 1001+291
is larger than the previous measurement of Du & Wang (2019) due to
variability.

We simply follow Sulentic et al. (2011) to define outliers in the Eigenvector
1 plane, namely, those AGNs with $\rm FWHM_{H\beta}>4000~{}km~{}s^{-1}$ and
${{\cal R}_{\rm Fe}}\,>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 ${\cal R}_{\rm Fe}$ and broad
H$\beta$ 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.

Figure 11: The main sequence of RM samples. The black points were compiled by
Du & Wang (2019). The blue horizontal and vertical lines corresponds to $\rm
FWHM_{H\beta}=4000~{}km~{}s^{-1}$ and ${{\cal R}_{\rm Fe}}=0.5$, respectively,
and the gray dotted line represent ${{\cal R}_{\rm Fe}}=1.0$. The upper right
corner enclosed by the blue solid lines is defined as the “outlier” region
according to Sulentic et al. (2011).

##  6 Conclusion

We present a reverberation mapping campaign of two luminous quasars at the
high-luminosity end of the $R_{\rm H\beta}-L_{5100}$ relation. Our main
results are as follows.

  * •

We measure the time delays of the broad emission lines with respect to the
5100 Å continuum. Using the ICCF method, the H$\beta$ and H$\gamma$ lags of PG
0923+201 in the rest frame are $108.2_{-12.3}^{+6.6}$ and
$121.0_{-10.2}^{+8.5}$ days, respectively, and the H$\beta$, H$\gamma$, Fe ii
lags of PG 1001+291 are $37.3_{-6.0}^{+6.9}$, $28.0_{-17.6}^{+17.9}$ and
$57.2_{-11.8}^{+10.5}$ days, respectively.

  * •

The velocity-resolved delays of the H$\beta$ line in PG 0923+201 show a
virialized motion, with shorter lags at line wings and longer lags at line
core. However, the lags in the redshifted bins are higher than those in the
blueshifted bins. This complicated structure may indicate the coexistence of
virialized motion and outflow in the BLR of PG 0923+201. The lag structure of
PG 1001+291 is relatively symmetric and also shows a virialized BLR.

  * •

Based on the H$\beta$ delays and FWHMs in the mean spectra, and assuming a
virial factor of $f_{\rm BLR}=1$, we estimate the black hole masses to be
$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. We obtain
consistent mass estimates using the H$\gamma$ line. The accretion rates of PG
0923+201 and PG 1001+291 are estimated to be $0.21_{-0.07}^{+0.06}\,L_{\rm
Edd}\,c^{-2}$ and $679_{-227}^{+259}\,L_{\rm Edd}\,c^{-2}$, respectively,
indicating that PG 0923+201 is accreting at a sub-Eddington rate whereas PG
1001+291 is a super-Eddington accretor.

  * •

The H$\beta$ time lag of PG 1001+291 falls 0.78 dex below the empirical
$R_{\rm H\beta}-L_{5100}$ relation of Bentz et al. (2013), confirming that
even for high-luminosity quasars, H$\beta$ time lags depend on both luminosity
and Eddington ratio, as previously found by Du et al. (2018b). This
strengthens the conclusion that the $R_{\rm H\beta}-L_{5100}$ relation at the
high-luminosity end needs to consider the influences of accretion rates. The
uncomfortably high single-epoch black hole masses estimated for AGN at large
redshifts may be significantly over-estimated if this effect has not been
taken into account.

We thank the referee for useful comments that improved the manuscript. This
work is supported by the National Key R&D Program of China (2016YFA0400701,
2016YFA0400702), by the National Science Foundation of China (NSFC-11721303,
11773029, 11833008, 11873048, 11922304, 11973029, 11991051, 11991052,
11991053, 11991054, 12003036, 12022301), by the Key Research Program of
Frontier Sciences of the Chinese Academy of Sciences (CAS; YZDJ-SSW-SLH007),
by the CAS Key Research Program (KJZDEW-M06), by the CAS International
Partnership Program (113111KYSB20200014), and by the Strategic Priority
Research Program of the CAS (XDB23000000, XDB23010400). We acknowledge the
support of the staff of the Lijiang 2.4 m telescope. Funding for the telescope
has been provided by the CAS and the People’s Government of Yunnan Province.
We also acknowledge the support of the staff of the CAHA 2.2 m telescope.
Y.-R.L. acknowledges financial support from the Youth Innovation Promotion
Association CAS. K.H. acknowledges support from STFC grant ST/R000824/1. Based
on observations obtained with the Samuel Oschin Telescope 48-inch and the
60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient
Facility project. ZTF is supported by the National Science Foundation and a
collaboration including Caltech, IPAC, the Weizmann Institute for Science, the
Oskar Klein Center at Stockholm University, the University of Maryland,
Deutsches Elektronen-Synchrotron and Humboldt University, Lawrence Livermore
National Laboratory, the TANGO Consortium of Taiwan, the University of
Wisconsin at Milwaukee, Trinity College Dublin, and Institut national de
physique nucléaire et de physique des particules. Operations are conducted by
COO, IPAC and University of Washington. This research has made use of the
NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion
Laboratory, California Institute of Technology, under contract with the
National Aeronautics and Space Administration.

Software: DASpec (https://github.com/PuDu-Astro/DASpec), PyCALI (Li et al.,
2014), MICA (Li et al., 2016), JAVELIN (Zu et al., 2011).

##  Appendix A Spectral Decompositions with Including the He ii Component

In this appendix, we present the fitting results of PG 0923+201 by adding a He
ii component. The SFS is similar to that described in Section 3.1 except for
adding a Gaussian to account for the He ii line (see Figure 12). In each
individual spectrum, the He ii is too weak and highly blended with the Fe ii.
Therefore, to reduce the degeneracy, the line width and shift of the He ii are
fixed to the best values obtained by fitting the rms spectrum (in the bottom
panel of Figure 3). The measured light curves and time lag analysis are shown
in Figure 13. The obtained time delays are summarized in Table 7.

Figure 12: A spectral decomposition of the mean spectrum for PG 0923+201 with including a He II $\lambda 4686$ component. The top and bottom panels show the details of the spectral decomposition and the residuals in the fitting window, respectively. The vertical dotted lines indicate the regions omitted from the fit. Table 7: Rest-frame time lag measurements by including a He ii component in spectral decompositions  Object | $\rm Line$ | $r_{\rm max}$ |  | $\tau_{\rm cent}$ | $\tau_{\rm MICA}$ | $\tau_{\rm JAV}$
---|---|---|---|---|---|---
| (day)
PG 0923+201 | ${\rm H\beta}$ | 0.92 |  | $110.1_{-10.5}^{+7.7}$ | $107.4_{-3.9}^{+4.5}$ | $107.8_{-3.4}^{+3.0}$
| ${\rm H\gamma}$ | 0.93 |  | $111.9_{-7.6}^{+10.2}$ | $107.6_{-4.4}^{+4.3}$ | $108.3_{-3.8}^{+3.9}$
Figure 13: Same as Figure 5, but for spectral decompositions including a He ii
component.

##  Appendix B ICCF Results in Each Velocity Bin

In Figure 14, we show the H$\beta$ light curves, ICCFs and CCCDs at different
velocity bins for PG 0923+201 and PG 1001+291.

Figure 14: The topmost panels show the continuum light curves and ACFs. The
rest panels show the H$\beta$ light curves (black points with error bars),
ICCFs (black lines) and CCCDs (blue histograms) at each velocity bin for
(left) PG 0923+201 and (right) PG 1001+291. The units of the continuum and
H$\beta$ light curves are $10^{-16}~{}\rm erg\,s^{-1}\,cm^{-2}$ Å-1 and
$10^{-14}~{}\rm erg\,s^{-1}\,cm^{-2}$, respectively. The velocity bin edges
are marked by pairs of numbers (in units of km s-1).

## References

  * Bahcall et al. (1972) Bahcall, J. N., Kozlovsky, B.-Z., & Salpeter, E. E. 1972, ApJ, 171, 467
  * Barth et al. (2015) Barth, A. J., Bennert, V. N., Canalizo, G., et al. 2015, ApJS, 217, 26
  * Barth et al. (2013) Barth, A. J., Pancoast, A., Bennert, V. N., et al. 2013, ApJ, 769, 128
  * Barth et al. (2011) Barth, A. J., Pancoast, A., Thorman, S. J., et al. 2011, ApJ, 743, L4
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002
  * Bentz et al. (2013) Bentz, M. C., Denney, K. D., Grier, C. J., et al. 2013, ApJ, 767, 149
  * Bentz et al. (2008) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2008, ApJ, 68
  * Bentz et al. (2009) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2009, ApJ, 705, 199
  * Bentz et al. (2010) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993
  * Blandford & McKee (1982) Blandford, R. D., & McKee, C. F. 1982, ApJ, 255, 419
  * Boroson & Green (1992) Boroson, T. A., & Green, R. F. 1992, ApJS, 80, 109
  * Brewer et al. (2011) Brewer, B. J., Páatay, L. B., & Csányi, G. 2011, Stat. Comput., 21, 649
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Czerny et al. (2019) Czerny, B., Wang, J.-M., Du, P., et al. 2019, ApJ, 870, 84
  * De Rosa et al. (2018) De Rosa, G., Fausnaugh, M. M., Grier, C. J., et al. 2018, ApJ, 866, 133
  * Denney et al. (2009) Denney, K. D., Peterson, B. M., Pogge, R. W., et al. 2009, ApJ, 704, L80
  * Denney et al. (2010) Denney, K. D., Peterson, B. M., Pogge, R. W., et al. 2010, ApJ, 721, 715
  * Du et al. (2014) Du, P., Hu, C., Lu, K.-X., et al. 2014, ApJ, 782, 45
  * Du et al. (2015) Du, P., Hu, C., Lu, K.-X., et al. 2015, ApJ, 806, 22
  * Du et al. (2016a) Du, P., Lu, K.-X., Hu, C., et al. 2016a, ApJ, 820, 27
  * Du et al. (2016b) Du, P., Lu, K.-X., Zhang, Z.-X., et al. 2016b, ApJ, 825, 126
  * Du et al. (2018a) Du, P., Brotherton, M. S., Wang, K., et al. 2018a, ApJ, 869, 142
  * Du & Wang (2019) Du, P. & Wang, J.-M. 2019, ApJ, 886, 42
  * Du et al. (2018b) Du, P., Zhang, Z.-X., Wang, K., et al. 2018b, ApJ, 856, 6
  * Edelson et al. (2019) Edelson, R., Gelbord, J., Cackett, E., et al. 2019, ApJ, 870, 123
  * Edelson et al. (2002) Edelson, R., Turner, T. J., Pounds, K. et al. 2002, ApJ, 568, 610
  * Fausnaugh et al. (2017) Fausnaugh, M. M., Grier, C. J., Bentz, M. C., et al. 2017, ApJ, 840, 97
  * Feng et al. (2021) Feng, H.-C., Hu, C., Li, S.-S., et al. 2021, ApJ, 909, 18
  * Fonseca Alvarez et al. (2020) Fonseca Alvarez, G., Trump, J. R., Homayouni, Y., et al. 2020, ApJ, 899, 73
  * Gaskell & Peterson (1987) Gaskell, C. M., & Peterson, B. M. 1987, ApJS, 65, 1
  * Gaskell & Sparke (1986) Gaskell, C. M. & Sparke, L. S. 1986, ApJ, 305, 175
  * Giveon et al. (1999) Giveon, U., Maoz, D., Kaspi, S., et al. 1999, MNRAS, 306, 637
  * Goodman & Weare (2010) Goodman, J. & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001
  * Grier et al. (2012) Grier, C. J., Peterson, B. M., Pogge, R. W., et al. 2012, ApJ, 755, 60
  * Grier et al. (2017) Grier, C. J., Trump, J. R., Shen, Y., et al. 2017, ApJ, 851, 21
  * Henden et al. (2012) Henden, A. A., Levine, S. E., Terrell, D., et al. 2012, Journal of the American Association of Variable Star Observers (JAAVSO), 40, 430
  * Ho & Kim (2014) Ho, L. C. & Kim, M. 2014, ApJ, 789, 17
  * Hook et al. (1994) Hook, I. M., McMahon, R. G., Boyle, B. J., et al. 1994, MNRAS, 268, 305
  * Hu et al. (2015) Hu, C., Du, P., Lu, K.-X., et al. 2015, ApJ, 804, 138
  * Hu et al. (2020a) Hu, C., Li, S.-S., Guo, W.-J., et al. 2020a, ApJ, 905, 75
  * Hu et al. (2021) Hu, C., Li, S.-S., Yang, S., et al. 2021, ApJS, 253, 20
  * Hu et al. (2020b) Hu, C., Li, Y.-R., Du, P., et al. 2020b, ApJ, 890, 71
  * Huang et al. (2019) Huang, Y.-K., Hu, C., Zhao, Y.-L., et al. 2019, ApJ, 876, 102
  * Kaspi et al. (2000) Kaspi, S., Smith, P. S., Netzer, H., et al. 2000,ApJ, 533, 631
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895
  * Kochanek et al. (2017) Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017, PASP, 129, 104502
  * Kollatschny et al. (2014) Kollatschny, W., Ulbrich, K., Zetzl, M., et al. 2014, A&A, 566, A106
  * Li et al. (2016) Li, Y.-R., Wang, J.-M., & Bai, J.-M. 2016, ApJ, 831, 206
  * Li et al. (2014) Li, Y.-R., Wang, J.-M., Hu, C., et al. 2014, ApJ, 786, L6
  * Lu et al. (2019) Lu, K.-X., Bai, J.-M., Zhang, Z.-X., et al. 2019, ApJ, 887, 135
  * Maoz et al. (1990) Maoz, D., Netzer, H., Leibowitz, E., et al. 1990, ApJ, 351, 75. doi:10.1086/168445
  * Martínez-Aldama et al. (2019) Martínez-Aldama, M. L., Czerny, B., Kawka, D., et al. 2019, ApJ, 883, 170
  * Marziani et al. (2001) Marziani, P., Sulentic, J. W., Zwitter, T., et al. 2001, ApJ, 558, 553
  * Marziani & Sulentic (2014) Marziani, P. & Sulentic, J. W. 2014, MNRAS, 442, 1211
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003
  * Mejía-Restrepo et al. (2018) Mejía-Restrepo, J. E., Lira, P., Netzer, H., et al. 2018, Nature Astronomy, 2, 63
  * Pei et al. (2017) Pei, L., Fausnaugh, M. M., Barth, A. J., et al. 2017, ApJ, 837, 131
  * Peterson (2014) Peterson, B. M. 2014, Space Sci. Rev., 183, 253
  * Peterson et al. (2002) Peterson, B. M., Berlind, P., Bertram, R., et al. 2002, ApJ, 581, 197
  * Peterson et al. (2004) Peterson, B. M., Ferrarese, L., Gilbert, K. M., et al. 2004, ApJ, 613, 682
  * Peterson (1993) Peterson, B. M. 1993, PASP, 105, 247
  * Peterson & Wandel (1999) Peterson, B. M. & Wandel, A. 1999, ApJ, 521, L95
  * Peterson & Wandel (2000) Peterson, B. M. & Wandel, A. 2000, ApJ, 540, L13
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Rodríguez-Pascual et al. (1997) Rodríguez-Pascual, P. M., Alloin, D., Clavel, J., et al. 1997, ApJS, 110, 9
  * Schlafly & Finkbeiner (2011) Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Shakura & Sunyeav (1973) Shakura, N. I. & Sunyaev, R. 1973, A&A, 24, 337
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48
  * Shen & Ho (2014) Shen, Y. & Ho, L. C. 2014, Nature, 513, 210
  * Shen et al. (2016) Shen, Y., Horne, K., Grier, C. J., et al. 2016, ApJ, 818, 30
  * Shen et al. (2011) Shen, Y., Richards, G. T., Strauss, M. A., et al. 2011, ApJS, 194, 45
  * Sulentic et al. (2000) Sulentic, J. W., Marziani, P., & Dultzin-Hacyan, D. 2000, ARA&A, 38, 521
  * Sulentic et al. (2011) Sulentic, J., Marziani, P., & Zamfir, S. 2011, Baltic Astronomy, 20, 427
  * van Groningen & Wanders (1992) van Groningen, E., & Wanders, I. 1992, PASP, 104, 700
  * Vanden Berk et al. (2004) Vanden Berk, D. E., Wilhite, B. C., Kron, R. G., et al. 2004, ApJ, 601, 692
  * Veilleux & Osterbrock (1987) Veilleux, S. & Osterbrock, D. E. 1987, ApJS, 63, 295
  * Wang et al. (2014a) Wang, J.-M., Du, P., Li, Y.-R., et al. 2014a, ApJ, 792, L13
  * Wang et al. (2013) Wang, J.-M., Du, P., Valls-Gabaud, D. & Netzer, H. 2013, Phys. Rev. Lett., 110, 081301
  * Wang et al. (2014b) Wang, J.-M., Qiu, J., Du, P., et al. 2014b, ApJ, 797, 65
  * Wang et al. (1999) Wang, J.-M., Szusszkiewicz, E., Zhou, Y.-Y. & Lu, F.-J. 1999b, ApJ, 522, 839
  * Wang & Zhou (1999) Wang, J.-M. & Zhou, Y.-Y. 1999a, ApJ, 516, 420
  * Yu et al. (2019) Yu, L.-M., Bian, W.-H., Wang, C., et al. 2019, MNRAS, 488, 1519
  * Zhang et al. (2019) Zhang, Z.-X., Du, P., Smith, P. S., et al. 2019, ApJ, 876, 49
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M., 2011, ApJ, 735, 80


Spiegelgasse 1
4051 Basel, Schweiz
Email<EMAIL_ADDRESS>