Broad-line Region of the Quasar PG 2130+099 from a Two-year Reverberation Mapping Campaign with High Cadence

The photometric images were obtained by CAFOS in direct imaging mode with a Johnson V filter. Typically, three exposures of 20 s were taken each individual night. The images were reduced following standard IRAF procedures. The flux of the object and the comparison star were measured through a circular aperture with a radius of 21, and differential magnitudes were obtained relative to 17 other stars within the field of 88 × 88.

Figure 1 shows the light curves of the differential magnitudes for the object (top) and the comparison star (bottom), respectively. The accurate of the photometry is ∼0.01 mag, estimated from the scatter in the magnitudes of the comparison star. It also shows that this comparison star is stable enough to be used for the flux calibration of the spectroscopy.

Zoom In Zoom Out Reset image size

The spectra of PG 2130+099 were taken using CAFOS with Grism G-200, which has a relatively high efficiency, without an order-blocking filter. The resulting spectra cover the observed-frame wavelength range of 4000–8500 Å, with a dispersion of 4.47 Å pixel⁻¹. Note that the second-order contamination occurs for wavelengths longer than ∼7000 Å, and has no effect on the measurements. A set of calibration frames was taken for each night, including bias frames, dome flats, and wavelength-calibration lamps of HgCd/He/Rb. For the object, two exposures of 600 s were taken during each individual night, and the typical signal-to-noise ratio (S/N) of the spectrum for a single exposure is ∼80 per pixel at rest frame 5100 Å. We took one or two spectrophotometric standards in the dusk and/or dawn of each night, if the weather and time allowed.

The spectra were reduced with IRAF following standard procedures: bias-removal, flat-fielding, wavelength calibration, and extraction to one-dimension. Extraction apertures were uniform and large (106) to minimize light loss, especially on nights with poor seeing. The flux calibration of the object used the sensitivity function determined from the comparison star, as described below. First, for several nights with good weather conditions, the normal IRAF procedure of flux calibration was performed for the comparison star using spectrophotometric standards taken on the same night. Second, the calibrated spectra for those nights were combined to generate a fiducial spectrum of the comparison star. Then, for each exposure, a Legendre polynomial was fitted by comparing the extracted spectrum of the comparison star in counts to the fiducial spectrum. This Legendre polynomial serves as the sensitivity function and incorporates all corrections including the atmosphere, slit loss, and instrumental sensitivity. Finally, the spectrum of the object for each exposure was flux-calibrated using the corresponding sensitivity function, and the two calibrated spectra for the same night were combined to obtain the individual-night spectrum for the following measurements and analysis.

The accuracy of the flux calibration using our comparison stars has been proven to be better than ∼3% (Du et al. 2018). In this specific case, it is better than 3% as estimated from the scatter of the [O iii] light curve (see Section 3.1.2 for details). Figure 2 shows the mean (top) and rms (bottom) spectra for the years 2017 (red) and 2018 (blue) separately. The [O iii] line almost vanishes in the rms spectrum, indicating that we have achieved good flux calibration. A prominent feature of the rms spectrum, compared to the mean spectrum, is the strong, broad He ii emission line, indicating the strong variation in this line. The mean spectra for the two years are almost identical, and the variations of the Hβ profiles are also quite small compared with those in Kaspi et al. (2000) and Grier et al. (2012). However, the rms spectrum of 2018 is dramatically different in its almost vanishing Balmer lines.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In 2018, the object was also monitored by two other telescopes: the Lijiang 2.4 m telescope and the Sutherland 1.9 m telescope. Spectra were taken on three nights at the Lijiang 2.4 m telescope at the Yunnan Observatory of the Chinese Academy of Sciences, using the Yunnan Faint Object Spectrograph and Camera with grism G14 and a long slit of 25. The slit was oriented to include the same comparison star as at CAHA, which is used for flux calibration. The spectra were extracted in an aperture of 25 × 85. The details of the spectroscopy and data reduction are given in Du et al. (2014). At the Sutherland station of the South African Astronomical Observatory, spectroscopic observations of PG 2130+099 were carried out with the 600 lines mm⁻¹ grating and a slit of 09, on 12 nights. The flux calibration was done by scaling the flux of the narrow [O iii] emission lines to be the same as the mean value in the CAHA spectra. The set of the observations and data reduction are the same as that described in Winkler & Paul (2017).

The epochs of the Lijiang spectra are rather few (only three nights), while the Sutherland spectra are calibrated by another method and the data quality is relatively scattered especially for the last four data points where technical problems caused a shift in the grating angle. So for the consistency of the data set, we used the data from CAHA only for the analysis below, but show the integrated light curves of the continuum and Hβ from Lijiang (green points) and Sutherland (blue points) in panels (a) and (c) of Figure 4.

Zoom In Zoom Out Reset image size

There are several methods to measure the light curves. The traditional method adopted in most RM studies (e.g., Kaspi et al. 2000; Peterson et al. 2004; Du et al. 2018) employs integration over a range of relevant wavelengths, which is simple and usually robust for single, strong emission lines such as Hβ. But for those highly blended emission lines, e.g., Fe ii emission which forms a pseudo-continuum, a spectral-fitting scheme (SFS) has proven to be necessary to measure the light curves (Bian et al. 2010; Barth et al. 2013; Hu et al. 2015). In some cases, especially for nearby Seyferts, spectral fitting can improve the measurements for even the Hβ emission line, by modeling and removing the contamination of strong host galaxy contribution which varies from night to night due to seeing and guiding variations (Hu et al. 2015, 2016). However, comparing with integration, the spectral fitting scheme requires higher quality of the data. For example, for reliable decomposition of the host starlight, the wavelength coverage of the spectrum should be wide, and the shape of the spectrum has to be especially well calibrated. These requirements can be difficult to achieve with the [O iii]-calibration mode of RM campaigns since only a constant scaling factor is applied for the whole wavelength range. Sometimes, integration and fitting are combined to use the advantages of both methods: subtracting the continuum by spectral fitting and then integrating the residual flux to measure the light curves of emission lines (Barth et al. 2013, 2015).

For fitting the spectra of PG 2130+099 here, besides an AGN power law, a host starlight template is needed to model the continuum shape accurately. However, the host starlight contribution is rather weak in this luminous PG object, and also the spectral resolution is low; our spectra do not show prominent absorption features of host starlight. Thus, the host starlight component cannot be modeled well and its apparent variation due to seeing and mis-centering cannot be removed by spectral fitting as precisely as for other sources (Hu et al. 2015, 2016). In fact, the uncertainty in fitting the host starlight could introduce extra systematic errors into the light-curve measurements. So in this work, we used both methods in this manner: integration for the measurements of the AGN continuum, Hβ, and He i, and the fitting method for the highly blended Fe ii and He ii features.

3.1.1. Integration Scheme

3.1.2. Spectral Fitting Scheme

We use the light curves of the continuum at 5100 Å, Hβ, and He i from integration, and He ii and Fe ii from spectral fitting, for the following time-series analysis. Table 1 lists all these light curves. The errors of the fluxes listed in this table are those generated from the errors of each pixel in the observed spectra, and are usually not large enough to interpret the scatter in the fluxes of successive nights. So an additional systematic error was calculated for each light curve using the same method as in Du et al. (2014), and all the analysis below was done taking account of this systematic error. We first calculate the variability amplitude ${F}_{\mathrm{var}}$, defined by Rodríguez-Pascual et al. (1997), for each light curve. The quantity represents the intrinsic variability over the errors in the measurements (including the additional systematic error). Its uncertainties are calculated as defined by Edelson et al. (2002). Table 2 lists the results. As mentioned in Section 3.1.2 above, the weak narrow Hβ component has not been subtracted from our measurements of broad Hβ. The flux of [O iii] λ5007 line is $1.21\pm 0.04\times {10}^{-13}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$, measured from the light curve of [O iii] in panel (d) of Figure 4. Assuming a typical value of 0.1 for the narrow Hβ/[O iii] intensity ratio in AGNs (e.g., Veilleux & Osterbrock 1987), the narrow Hβ component contributes only ∼2% fluxes of the Hβ measured here. Thus, the dilution of Hβ ${F}_{\mathrm{var}}$ by the uncorrected narrow component is negligible. Note that He ii has an ${F}_{\mathrm{var}}$ much larger than the other emission lines, and even the continuum. This behavior is also shown by the strong He ii emission feature in the rms spectrum (see the bottom panel of Figure 2), and common for many objects in previous RM investigations (Barth et al. 2015).

Table 1.  Light Curves of 5100 Å Continuum and Several Lines

JD-2457,900	${F}_{5100}$	${F}_{{\rm{H}}\beta }$	${F}_{\mathrm{He}{\rm{I}}}$	${F}_{\mathrm{He}{\rm{II}}}$	${F}_{\mathrm{Fe}}$
27.637	5.314 ± 0.017	5.927 ± 0.010	0.590 ± 0.006	0.610 ± 0.015	4.636 ± 0.020
31.593	5.298 ± 0.037	5.993 ± 0.012	0.587 ± 0.007	0.459 ± 0.018	4.569 ± 0.025

Note. The 5100 Å continuum flux is in units of ${10}^{-15}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathring{\rm A} }^{-1}$, and all lines are in units of ${10}^{-13}\,\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 2.  Reverberations of Several Broad Lines from the Present Campaign

Line	Fvar (%)		Lag (days)		FWHM ($\mathrm{km}\,{{\rm{s}}}^{-1}$)		Virial Product $(\times {10}^{7}{M}_{\odot })$		CCF Broadening (days)
	2017	2018	2017	2018	2017	2018	2017	2018	2017	2018
He ii	32.7 ± 3.6	29.9 ± 3.0	$-{1.4}_{-0.9}^{+6.8}$	$-{2.9}_{-0.9}^{+1.5}$	6049 ± 187	9589 ± 774	⋯	⋯	5	−7
He i	8.2 ± 0.9	4.2 ± 0.6	${18.2}_{-3.1}^{+7.3}$	${31.1}_{-5.8}^{+2.9}$	2547 ± 18	2497 ± 14	${2.30}_{-0.39}^{+0.93}$	${3.78}_{-0.70}^{+0.35}$	17	7
Hβ	6.1 ± 0.6	2.3 ± 0.3	${22.6}_{-3.6}^{+2.7}$	${27.8}_{-2.9}^{+2.9}$	2101 ± 100	2072 ± 107	${1.95}_{-0.36}^{+0.30}$	${2.33}_{-0.34}^{+0.34}$	31	12
Fe ii	4.8 ± 0.5	3.3 ± 0.5	${35.3}_{-9.9}^{+8.2}$	${23.1}_{-5.6}^{+3.4}$	2047 ± 94	1944 ± 95	${2.88}_{-0.85}^{+0.72}$	${1.70}_{-0.45}^{+0.30}$	48	16

Note. The virial product is defined by Equation (1). See the text for details of the measurements of FWHM for each line. It should be pointed out that ${F}_{\mathrm{var}}$ has different meanings from the calibration uncertainties.

Download table as:  ASCIITypeset image

We measured the reverberation lags between the variations of the continuum (${F}_{5100}$) and emission lines (${F}_{\mathrm{He}{\rm{II}}}$, ${F}_{\mathrm{He}{\rm{I}},\mathrm{intg}}$, ${F}_{{\rm{H}}\beta ,\mathrm{intg}}$, and ${F}_{\mathrm{Fe}}$), using the standard interpolation cross-correlation function (CCF) method (Gaskell & Sparke 1986; Gaskell & Peterson 1987; White & Peterson 1994). The centroid of the CCF, using only the part above 80% of its peak value (${r}_{\max }$), is adopted as the time lag (Koratkar & Gaskell 1991; Peterson et al. 2004). The error of the time lag is given by the 15.87% and 84.13% quantiles of the cross-correlation centroid distribution (CCCD) generated by 5000 Monte Carlo realizations, following Maoz & Netzer (1989) and Peterson et al. (1998b). In each realization, a subset of data points is selected randomly, and their fluxes are also changed by random Gaussian deviations according to the errors. The top panel of Figure 5 shows the autocorrelation function (ACF) of the ${F}_{5100}$ light curve (in black), and also the CCFs for the emission lines (He ii in cyan, He i in red, Hβ in green, and Fe ii in blue) with respect to ${F}_{5100}$. The bottom panel shows the corresponding CCCDs for the CCFs. The measured time lags are listed in Table 2.

Zoom In Zoom Out Reset image size

For Hβ and Fe ii, which are included in the fitting and allowed to vary, the means of their velocity widths and shifts obtained from the best fits to individual-night spectra are adopted to be the measurements of their line profiles. The standard deviations are used as the errors. For He ii, the velocity width and shift are measured from the rms spectrum. For He i, we perform an additional fit to a narrow band of the spectrum around the emission line including only a straight line as the continuum and a single Gaussian for He i. Their errors are estimated by Monte Carlo simulation: several realizations of mean and rms spectra were generated by bootstrap sample selection, and the same spectral fitting was performed and the velocity widths were measured. Then the standard deviations of the velocity widths measured from the realizations are used as the errors.

The resolution is estimated by comparing the width of [O iii] in our spectra (FWHM = 1084 km ${{\rm{s}}}^{-1}$) to those from previous high-spectral-resolution observations. Whittle (1992) obtained an [O iii] FWHM measurement of 350 km ${{\rm{s}}}^{-1}$ (see also Grier et al. 2012), yielding a broadening of 1020 km ${{\rm{s}}}^{-1}$. We also fit the spectrum of this object from Boroson & Green (1992), obtained an [O iii] FWHM of 400 km ${{\rm{s}}}^{-1}$ after correcting for their instrument broadening. The resultant broadening is similar, and we adopted an FWHM = 1000 km ${{\rm{s}}}^{-1}$ as the instrument broadening of our spectra. The measurements of the widths of each emission line after instrument broadening correction are listed in Tables 2 and 3.

Table 3.  Hβ Line Width Information

	Mean Spectra			Rms
	σline			FWHM			σline
Year	2017	2018		2017	2018		2017	2018
Width	1437 ± 64	1438 ± 41		2447 ± 73	2054 ± 557		1272 ± 97	874 ± 237
${\hat{M}}_{\bullet }({10}^{7}{M}_{\odot })$	${0.91}_{-0.17}^{+0.14}$	${1.12}_{-0.13}^{+0.13}$		${2.64}_{-0.45}^{+0.35}$	${2.29}_{-1.26}^{+1.26}$		${0.71}_{-0.16}^{+0.14}$	${0.41}_{-0.23}^{+0.23}$

Note. The width is in units of km ${{\rm{s}}}^{-1}$ and after instrumental broadening correction. Here the values were measured from the rms spectra shown in Figure 2, without the deconvolution in Section 4.3. The blue peak around −2500 km ${{\rm{s}}}^{-1}$ in the rms spectrum after deconvolution in Figure 6 is not high enough to contribute the FWHM listed here.

Download table as:  ASCIITypeset image

In order to investigate the geometry and kinematics of the BLR in PG 2130+099, we calculated the velocity-resolved time lags of the Hβ emission line profile. First, the influence of the varying line-broadening functions (ψ) and wavelength-calibration inaccuracy for different nights, which are caused by the changing seeing and misalignment of the slit, must be removed. Using comparison stars, the line-broadening function can be obtained by fitting the spectra of the stars using a stellar template convolved by ψ. The wavelength-calibration inaccuracy is also taken into account as the velocity shift of ψ. We adopted the Richardson–Lucy deconvolution algorithm, which is demonstrated to be an efficient way to recover the signal blurred by a known response kernel, to correct the ψ function. Details of these mathematical manipulations are described by Du et al. (2016a). After the correction, the emission line is divided into several bins, each of which has the same flux in the rms spectrum. Then, the light curves in the bins and their corresponding time lags relative to the continuum light curve are measured as previously described.

Figure 6 shows the velocity-resolved time delays and the rms spectrum (after the line-broadening correction) obtained for 2017 (left) and 2018 (right), respectively. The figure indicates that the BLR in 2017 has a geometry of virialized motion (e.g., see plots in Welsh & Horne 1991; Bentz et al. 2009), while in 2018 it more likely has an infall structure. However, it should be pointed out that the bluest point in the right panel of Figure 6 has quite large error bars because of the weak response in the second year. Whether or not the BLR kinematics are consistent with inflow remains open, but it is certain that the variable BLR has changed between the two years. For the upcoming discussion, we refer to the BLR as an inflow in the second year. The exact description of the geometry should be given by the application of the maximum-entropy method (MEM) developed by Horne et al. (2004), which is out of the scope of this paper.

Zoom In Zoom Out Reset image size

We note that the lags in both wings could be contaminated by the continuum subtraction. In order to test the contamination, we show the velocity-resolved delays with equal fluxes in the mean spectra in Appendix A. Light curves of each velocity bin are shown in Appendix B as well as CCF analysis. Comparing them, we conclude that the influence is quite weak.

The data quality of other emission lines is not good enough to determine the velocity-resolved time lags.

Following the procedure in Li et al. (2016),¹³ we used a Gaussian to parameterize the transfer function of each emission line. The amplitude, width, and central values of the Gaussian are free parameters and determined by a Markov Chain Monte Carlo (MCMC) method. The continuum variations are described by a damped random walk model. As such, the covariance functions between continuum and emission lines and between lines can be expressed analytically. This allowed us to employ the well-established framework given by Rybicki & Press (1992) to calculate the Bayesian posterior probabilities and infer the best values for free parameters.

Fitting results for the two years are shown by Figure 7 and listed in Table 4, separately. Lags obtained by the transfer functions are consistent with those of the CCF method.

Zoom In Zoom Out Reset image size

The bulge velocity dispersion of PG 2130+099 is σ_* = 163 ± 19 km ${{\rm{s}}}^{-1}$ from near-infrared spectra (Grier et al. 2013a), but the host galaxy has a pseudo-bulge (n ≈ 0.45) and a disturbed disk (Kim et al. 2017). Using the same slope of the ${M}_{\bullet }\mbox{--}{\sigma }_{* }$ relation (Kormendy & Ho 2013) but half of its zero-point (Ho & Kim 2015), we expect ${M}_{\bullet }\approx 6.5\times {10}^{7}{M}_{\odot }$ in PG 2130+099. This estimate should be regarded as an upper limit. Indeed, Grier et al. (2017) made use of an MCMC method to model the BLR and to determine the black hole mass through Hβ light curves, and found ${\mathrm{log}}_{10}({M}_{\bullet }/{M}_{\odot })={6.87}_{-0.23}^{+0.24}$, which is significantly lower than that from the ${M}_{\bullet }\mbox{--}\sigma$ relation. Comparing our velocity-resolved time lags in 2017 as shown by Figure 6 with the MCMC model (Grier et al. 2017), we find the model matches our measurements, leading to support of ${M}_{\bullet }={10}^{6.87}{M}_{\odot }$ as a conservative mass estimation.

Using the presently detected Hβ lags, we determine the virial product

$$\begin{eqnarray}&&{\hat{M}}_{\bullet }=\displaystyle \frac{c{\tau }_{{\rm{H}}\beta }\times {V}^{2}}{G},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, V is either full width at half maximum (V_FWHM) or velocity dispersion (σ_line) of the Hβ profile in the mean spectrum or rms. We list ${\hat{M}}_{\bullet }$ in Tables 2 and 3, yielding the black hole mass if given the virial factor (${f}_{\mathrm{BLR}}$). The ${f}_{\mathrm{BLR}}$ can in principle be calibrated by the ${M}_{\bullet }\mbox{--}{\sigma }_{* }$ relation (e.g., Onken et al. 2004; Woo et al. 2015; Batiste et al. 2017), or the black hole masses from accretion disk models (Lü 2008; Mejía-Restrepo et al. 2017), but it tends to be smaller for AGNs with pseudo-bulges, such as ${f}_{\mathrm{BLR}}=0.5$ (Ho & Kim 2014). If we take ${f}_{\mathrm{BLR}}=0.5$, we have ${M}_{\bullet }={0.97}_{-0.18}^{+0.15}\times {10}^{7}\,{M}_{\odot }$ from the virial product by using the mean FWHM of Hβ in 2017 (the BLR may significantly deviate from virialized state in the second year), agreeing with the MCMC result of black hole mass in Grier et al. (2017).

In light of the standard model of accretion disks (Shakura & Sunyaev 1973), the dimensionless accretion rates defined by $\dot{{\mathscr{M}}}={\dot{M}}_{\bullet }/{L}_{\mathrm{Edd}}{c}^{-2}$ can be estimated as (Du et al. 2016b)

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

if given the optical luminosity and the black hole mass, where ${\dot{M}}_{\bullet }$ is the accretion rate, and L_Edd is the Eddington luminosity, ${{\ell }}_{44}={L}_{5100}/{10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, ${M}_{7}={M}_{\bullet }/{10}^{7}{M}_{\odot }$, and i is the inclination of the disks (we take cos i = 0.75¹⁴ ). The mean flux at 5100 Å in the first year is ${\bar{F}}_{5100}=(5.12\pm 0.25)\,\times {10}^{-15}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathring{\rm A} }^{-1}$ corresponding to a mean luminosity $\lambda {L}_{\lambda }=(2.50\pm 0.12)\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ in the present epoch. With ${M}_{\bullet }$ and ${{\ell }}_{44}=2.50$, we obtain $\dot{{\mathscr{M}}}={10}^{2.1\pm 0.5}$, implying a super-Eddington accretor.

The onion structure of the BLR is quite well understood as a sequence of ionization energy of ions (Collin-Souffrin & Lasota 1988) and FWHM of lines (Peterson & Wandel 1999). The ions relevant to this study have ionization energies of ${\epsilon }_{\mathrm{He}{\rm{II}}}=54.4\,\mathrm{eV}$, ${\epsilon }_{\mathrm{He}{\rm{I}}}=24.58\,\mathrm{eV}$, ${\epsilon }_{{\rm{H}}\beta }=13.6$ eV and ${\epsilon }_{\mathrm{Fe}\mathrm{II}}=7.87\,\mathrm{eV}$. From the first year, we have ${R}_{\mathrm{Fe}{\rm{II}}}\,\gt {R}_{{\rm{H}}\beta }\gt {R}_{\mathrm{He}{\rm{I}}}\gt {R}_{\mathrm{He}{\rm{II}}}$, following the onion structure as shown by Figure 8.

Zoom In Zoom Out Reset image size

In the observations of the second year, the Hβ-emitting region do not show significant variations; however, He i and Fe ii regions show significant changes, in particular, ${R}_{\mathrm{Fe}\mathrm{II}}\lesssim {R}_{{\rm{H}}\beta }\lesssim {R}_{\mathrm{He}{\rm{I}}}$ disrupting the onion structure with a timescale much less than the dynamical time as shown by Equation (3). The first unresolved question would be to explain the mechanism that generates such a change in the onion structure.

The velocity-resolved delays of the broad Hβ line in the two years shown in Figure 6 are different. The responding emission line region shows a change from a virialized motion to an inflow on a timescale of less than one year. The profile and the intensity of the Hβ line in the rms spectra dramatically change. It is a smooth core-dominated profile in the first year, but multiple-peaked and wing-dominated in the second. On the other hand, the widths of the broad lines in the mean spectrum remain almost the same as in the first year, as well as the virial products derived. This implies that the BLR may still be actually virialized in the second year. Such an inconsistency is worth further investigation: is the velocity-resolved delays in the second year a unique signature of infall as in those simple models? If so, why do the majority of Hβ-emitting clouds have virialized motions in the second year as shown by the mean spectrum, while the reverberation part has a kinematics of inflow indicated by the velocity-resolved delays? Furthermore, what is the connection between the inflow and the virialized part of the BLR? This is the second puzzle, which motivates us to continue monitoring PG2130+099.

The variability amplitude of 5100 Å flux (F_var = 6.7%) in the second year is larger than that in the first year (F_var = 4.7%); however, Hβ is less variable in the second year (F_var = 2.3%) than in the first year (F_var = 6.1%). This means a weaker response of Hβ to the varying continuum in the second year. That the virialized part of the BLR did not respond in the second year may be consistent with the weak response, but it is hard to understand why. The transfer functions are also very different in the two years. Though the mean radius of Hβ-emitting region does not change much, the width (τ/Δ τ) changes by a factor of 2 (see Table 4 and Figure 7). On the contrary, He i and Fe ii mean radii change significantly, but their relative widths (τ/Δτ) stay almost unchanged. This is the third puzzle.

Photoionization calculations and light-curve simulations have been performed in the literature (e.g., Korista & Goad 2004; Goad & Korista 2014; Lawther et al. 2018) to investigate how the emission line responsivity depends on factors, including the continuum state, driving continuum variability, BLR geometry, and the duration/cadence of the campaign. Korista & Goad (2004) calculated photoionization models suitable for NGC 5548, and found that the local line responsivity could be two times lower in the high continuum state than that in the low state (their Figure 3) with a factor of ∼8 times change in the continuum flux. In this case, the Hβ responsivities, estimated using the ${F}_{\mathrm{var}}$ ratio of the line to the continuum as suggested in Goad & Korista (2014), are ∼1.3 and 0.34 for the two years, respectively. But the continuum flux F₅₁₀₀ is only slightly higher (<5% on average, see Figure 4 panel (a)) in the second year. Note that the ionizing UV continuum may has larger variability between the two years than F₅₁₀₀. Its amplitude and impact on the responsivity need to be confirmed by further observations and calculations. Goad & Korista (2014) show that continuum variations faster than the maximum time lag for an extended BLR reduce the measured line responsivity and also the time lag in the meantime (their Figure 9). The continuum variability timescale in the second year seems to be somewhat faster than that in the first year by comparing the widths of the two ACFs in Figure 5, but the measured Hβ time lag does not become shorter accordingly as expected if the variability timescale is the reason for the weaker response. Thus, factors other than the continuum flux state and variability timescale are needed to interpret the much weaker responsivity in the second year. Detailed calculations and simulations out of the scope of this paper are necessary to answer this question.

As mentioned previously, three campaigns have been done for PG 2130+099, including Kaspi et al. (2000) and Grier et al. (2008, 2012) before the present work. Since the Kaspi et al. (2000) campaign has poor sampling and the Grier et al. (2008) campaign is too short in duration, we only compare our results with Grier et al. (2012), which performed a cadence of roughly one day. Grier et al. (2012) obtained ${\tau }_{{\rm{H}}\beta }=9.7\pm 1.3$ days, which is shorter than the value measured in this campaign by a factor of larger than two. Considering that our campaign is 7 yr later than that of Grier et al. (2012), the difference could be real and due to changes of the BLR. However, Bentz et al. (2013) argued that the Grier et al. (2012) campaign missed some key observation dates and PG 2130+099 should have a lag of 31 ± 4 days from their reanalysis. This corrected time lag is in agreement with the present results.