The Seoul National University AGN Monitoring Project. III. Hβ Lag Measurements of 32 Luminous Active Galactic Nuclei and the High-luminosity End of the Size–Luminosity Relation

We present the main results from a long-term reverberation mapping campaign carried out for the Seoul National University AGN Monitoring Project (SAMP). High-quality data were obtained during 2015–2021 for 32 luminous active galactic nuclei (AGNs; i.e., continuum luminosity in the range of 1044–46 erg s−1) at a regular cadence, of 20–30 days for spectroscopy and 3–5 days for photometry. We obtain time lag measurements between the variability in the Hβ emission and the continuum for 32 AGNs; 25 of those have the best lag measurements based on our quality assessment, examining correlation strength and the posterior lag distribution. Our study significantly increases the current sample of reverberation-mapped AGNs, particularly at the moderate-to-high-luminosity end. Combining our results with literature measurements, we derive an Hβ broadline region size–luminosity relation with a shallower slope than reported in the literature. For a given luminosity, most of our measured lags are shorter than the expectations, implying that single-epoch black hole mass estimators based on previous calibrations could suffer large systematic uncertainties.


Introduction
Black hole mass (M BH ) is a key parameter in understanding the physics of active galactic nuclei (AGNs) and the connection of black hole growth with galaxy evolution.M BH can be determined based on spatially resolved data by measuring the kinematics of stars, gas, and masers near the sphere of influence of supermassive black holes (e.g., Barth et al. 2001;Marconi et al. 2003;Davies et al. 2006;Gültekin et al. 2009;Scharwächter et al. 2013;den Brok et al. 2015;Boizelle et al. 2021;Kabasares et al. 2022) or by imaging black hole shadows along with a theoretical approach (Event Horizon Telescope Collaboration et al. 2019Collaboration et al. , 2022)).
However, these methods are limited to relatively nearby objects due to the limited spatial resolution of current facilities.
In contrast, reverberation mapping (RM; Blandford & McKee 1982;Peterson 1993) based on the variability of AGNs can be applied to mass-accreting black holes beyond the local Universe.Currently, the RM technique and related indirect mass estimators are the primary methods for determining M BH over a large range of cosmic time.The main idea of RM is to measure the time delay (τ) between the flux variations of the continuum and broad emission lines, which represents the light travel time from the central photoionizing source to the photoionized gas, providing the size (or radius) of the broadline region (BLR).Based on the virial assumption that the dynamics of gas in the BLR is governed by the gravitational potential of the central black hole, M BH is determined by combining the measured size (R BLR ) with the velocity measure (V ) from the width of broad emission lines as where c is the speed of light, G is the gravitational constant, and f is a factor representing the unknown geometry of the BLR.While f can be different for each AGN, an average f factor is calibrated based on the black hole mass-stellar velocity dispersion relation of the local galaxies (e.g., Onken et al. 2004;Woo et al. 2010Woo et al. , 2015;;Park et al. 2012a).Note that the f factor is the main source of the uncertainty of M BH , up to 0.4 dex in the case of Hβ-reverberation-based M BH (Park et al. 2012a;Woo et al. 2015).The f factor has been constrained for a small number of individual AGNs based on the dynamical modeling of the BLR combined with the velocity-resolved RM data (e.g., Pancoast et al. 2011Pancoast et al. , 2014aPancoast et al. , 2014b;;Li et al. 2013Li et al. , 2018;;Williams et al. 2018;Villafaña et al. 2022).
Early studies of RM reported a correlation between the measured Hβ BLR size and the monochromatic continuum luminosity at 5100 Å (L 5100 ; Wandel et al. 1999;Kaspi et al. 2000), opening an indirect way of estimating BLR size and M BH from single-epoch spectra, since monitoring data for RM are not required (e.g., Woo & Urry 2002;Vestergaard & Peterson 2006;Shen et al. 2011;Shen & Liu 2012;Dalla Bontà et al. 2020).While the best-fit slope was initially reported as 0.7 (Kaspi et al. 2000), following studies based on Hubble Space Telescope (HST) images calibrated the slope as ∼0.5, as expected from photoionization, after correcting for the host galaxy contribution to the observed L 5100 , particularly for lowluminosity AGNs with a relatively strong stellar continuum (e.g., Bentz et al. 2009aBentz et al. , 2013)).
It is important to investigate the photoionization and the BLR size-luminosity relation for a general population of AGNs over a large dynamic range of M BH and AGN luminosity.In the past, however, RM studies were limited to low-to-moderateluminosity AGNs, due to the observational challenges.The main difficulty is that long-term spectroscopic monitoring with a good cadence is required to obtain an accurate time lag measurement between AGN continuum and emission-line flux variations.Over the last decades, a number of intensive programs have been dedicated to RM studies, dramatically increasing the sample size and the dynamic range of the reverberation-mapped AGNs (Barth et al. 2015;Du et al. 2015;Grier et al. 2017;Zhang et al. 2019;U et al. 2022;Malik et al. 2023).However, it is still important to extend the RM study to more luminous AGNs, particularly in the high-luminosity regime (i.e., L 5100 = 10 45-46 erg s −1 ), which is the representative luminosity of high-z AGNs.For example, an AGN with L 5100 = 10 46 at z = 1 is expected to have an Hβ time lag of 500-600 days in the observed frame, which then has to be determined based on a monitoring campaign of 5-10 yr.Such a long time line explains why there is a relative lack of very luminous AGNs among the reverberation-mapped AGNs (see Figure 1).
Currently, the Hβ time lag has been reported for more than 200 AGNs (e.g., Wandel et al. 1999;Kaspi et al. 2000;Peterson et al. 2004;Bentz et al. 2009bBentz et al. , 2013;;Barth et al. 2015;Grier et al. 2017;Park et al. 2017;Du & Wang 2019;Martínez-Aldama et al. 2019;Zhang et al. 2019;Dalla Bontà et al. 2020;Hu et al. 2021;Li et al. 2021;U et al. 2022;Malik et al. 2023).These AGNs show a larger scatter in the Hβ BLR size-luminosity relation compared to the previously reported relations (e.g., Bentz et al. 2013).While it is possible that the intrinsic scatter may not be larger than previously thought, a consistent study of cross-correlations with uniform measurements of Hβ lag and uncertainty is required to properly constrain the intrinsic scatter as well as the slope of the sizeluminosity relation.Note that various studies performed by different groups have adopted different criteria to select reliable lag measurements and inconsistent methods to derive the uncertainty of the lag.
In particular, AGNs with a super-Eddington ratio seem to deviate from the size-luminosity relation according to the results from a recent project, the Super-Eddington Accreting Massive Black Holes (SEAMBHs; Du et al. 2014;Wang et al. 2014;Zhang et al. 2019;Hu et al. 2021;Li et al. 2021).It is claimed that the deviation from the size-luminosity relation correlates with the accretion rate or the flux ratios between Fe II and Hβ (R Fe ) or between [O III] and Hβ lines.This systematic trend is crucial to verify, since it will introduce strong bias in M BH determination for high-Eddington-ratio AGNs from single-epoch spectra (see the discussion by Li et al. 2021).
To extend the RM study to moderate-to-high-luminosity AGNs and investigate the Hβ BLR size-luminosity relation at the high-luminosity end, we performed an intensive long-term campaign, the Seoul National University AGN Monitoring Project (SAMP).Using 100 AGNs with moderate to high luminosity (L 5100 > 10 44.0 erg s −1 ) out to z < 0.4, we started test photometry in 2015.Then we carried out photometric and spectroscopic monitoring with the Lick 3 m, MDM 2.4 m, and other 1 m class telescopes until the middle of 2021.The project strategy and sample selection were presented by (Woo et al. 2019;hereafter, Paper I), and initial measurements of the Hβ lag for two targets based on the first 3 yr data were reported by Rakshit et al. (2019;hereafter Paper II).In this work, we present the final sample of 32 AGNs from the 6 yr spectroscopic campaign and the Hβ lag analysis.The measurements of Hα lag are presented by Cho et al. (2023;Paper IV), and the photometry monitoring results of 72 AGNs will be presented  (Bentz et al. 2009b;U et al. 2022), the monitoring campaign of PG quasars by Kaspi et al. (2000), SEAMBHs (Du et al. 2014(Du et al. , 2015(Du et al. , 2016(Du et al. , 2018a;;Wang et al. 2014;Zhang et al. 2019;Hu et al. 2021;Li et al. 2021), OzDES (Malik et al. 2023), SDSS-RM (Grier et al. 2017), and other AGNs in the collection by Dalla Bontà et al. (2020).by D. Son et al. (2024, in preparation).We briefly describe the sample selection in Section 2 and the data analysis in Section 3. Results and discussion are presented in Sections 4 and 5, respectively, followed by a discussion in Section 5 and a summary in Section 6.Throughout this paper, we use the ΛCDM cosmology, with H 0 = 72.0km s −1 Mpc −1 and Ω m = 0.3.

Sample
We selected the best available type 1 AGNs for a multiyear monitoring program for our facilities, by considering the expected lag, observability, and feasibility of Hβ lag measurements, along with the simulation of light curves and spectral decomposition.The details of the project strategy and sample selection of SAMP were presented in Paper I, and here we briefly describe the sample.
We selected relatively high-luminosity AGNs using the MILLIQUAS catalog (Milliquas v4.5 update ;Flesch 2015), in order to test the high end of the Hβ BLR size-luminosity relation.As summarized in Paper I, we initially selected 100 AGNs with V-band magnitude <17.0 at z < 0.5.The Hβ lag is expected to range from ∼40 to ∼250 days in the observed frame, which is estimated based on the monochromatic continuum luminosity at 5100 Å (L 5100 ) and the sizeluminosity relation of Bentz et al. (2013).In this process, we used Sloan Digital Sky Survey (SDSS) spectra or the spectra provided by Boroson & Green (1992) to measure λL λ (5100Å).B-band photometry is used instead if there is no available spectrum.We identified 48 AGNs with an expected lag longer than 70 days in the initial sample as the first priority targets (i.e., SAMP IDs starting with Pr1 in Table 1) and 37 AGNs with an expected lag shorter than 70 days as the second priority targets (i.e., SAMP IDs starting with Pr2).We also included 15 PG QSOs (Boroson & Green 1992) for filling up the seasonal gaps and increasing the sample size (i.e., SAMP IDs starting  with P).These PG QSOs are also medium-to-high-luminosity objects with an expected lag of ∼50-300 days.
Using the initial sample of 100 AGNs, we performed the variability test of continuum and Hβ line emission based on photometry and spectroscopy during the first few years, and excluded weakly varying objects.Based on these variability check processes, we continued 6 yr monitoring for a final sample of 32 AGNs.This final AGNs are out to z ∼ 0.4 with luminosity L 5100 > 10 44 erg s −1 , as presented in Figure 1.We emphasize that the SAMP final sample covers relatively higher luminosity ranges compared to previous Hβ monitoring campaigns (e.g., Kaspi et al. 2000;Bentz et al. 2009b;Du et al. 2014Du et al. , 2015Du et al. , 2016;;Wang et al. 2014;Barth et al. 2015;Grier et al. 2017;U et al. 2022;Malik et al. 2023).The properties of the sample are summarized in Table 1.

Photometry
We performed photometric monitoring using several telescopes: the MDM 1.3 and 2.4 m telescopes at Kitt Peak, Tucson, Arizona, USA; the Lemmonsan Optical Astronomy Observatory 1 m telescope located on Mount Lemmon, Tucson, Arizona, USA; the Lick Observatory 1 m telescope located at Mount Hamilton, California, USA; the Las Cumbres Observatory Global Telescope (LCOGT) network; and the Deokheung Optical Astronomy Observatory (DOAO) 1 m telescope.We used the B and V, or V and R, band filters to monitor continuum variability, depending on the redshift of each object (see Table 1).
Details of the photometric data reduction and variability analysis, including the final photometric light curves, will be presented in a forthcoming paper (D. Son et al. 2024, in preparation).Here we briefly describe the basic information for completeness.We followed the standard reduction procedure for bias subtraction and flat-fielding using the IRAF package (Tody 1986(Tody , 1993)).We used the LA-Cosmic 15 task (van Dokkum 2001) to remove cosmic rays and adopted   (Bertin et al. 2002), by matching the positions of all stars in the field of view (FOV) of each image.We performed aperture photometry using the SExtractor18 software (Bertin & Arnouts 1996), with an initial aperture size 3 times the seeing size.By generating the magnitude growth curve as a function of the aperture size, we tested whether the magnitude based on the seeing was fainter than the brightest magnitude in the growth curve, and enlarged the aperture size accordingly in order to avoid any aperture loss.A small fraction of images with bad quality due to a full moon, gust wind, or thick clouds were excluded from the final photometric light curves, based on the visual inspection of each image.We performed flux calibration using nonvariable stars in the FOV of each AGN, using the star catalogs of SDSS and APASS. 19inally, we obtained the photometric light curves of the sample in each band at each telescope.We found systematic offsets among the light curves obtained with different telescopes, presumably caused by various weather conditions and differences in the filter properties, as often reported by previous studies (e.g., Peterson et al. 1995;Pancoast et al. 2019).Thus, it was critical to intercalibrate these light curves in the merging process.We performed the intercalibration by adopting the Python software PyCALI20 (Li et al. 2014), which used a damped random walk (DRW) model to describe the AGN variability and determine the best scaling factors based on Bayesian statistics.For each AGN, we adopted the light curve from the MDM 1.3 m telescope as a reference and aligned all the other light curves.Note that the MDM 1.3 m light curves have the largest number of epochs, which are also most evenly distributed over the monitoring time baseline.Systematic uncertainties were added to each telescope's light curve during the intercalibration process.
The visual inspection of the merged light curves showed that for most targets, the DOAO light curves were still relatively scattered and deviated from the general trend of the light curves obtained at other telescopes, presumably due to flux calibration issues, i.e., high humidity and quickly changing seeing conditions at DOAO.In the case of the light curves obtained at the Lick 1 m and LCOGT telescopes, we sometimes found a similar problem for some objects.Under these circumstances, we decided to exclude the light curves from Lick 1 m and LCOGT as well.
We utilized the B-band light curves as the continuum light curves for all objects except for two higher-redshift objects at z > 0.3, namely J1105+671 and PG 2251+113, for which we instead adopted the V-band light curve.For seven AGNs, we found that the g-band light curves from the Zwicky Transient Facility (ZTF) were available, and combined them with the SAMP continuum light curves in order to improve the temporal coverage and cadence.Note that ZTF is a time-domain survey, which started in 2018 and overlapped with the SAMP monitoring baseline.We used ZTF Data Release 8 to include the photometric data from 2018 March to 2021 September.Following the previous studies based on ZTF data (e.g., Sánchez-Sáez et al. 2021), we cleaned the ZTF light curves by requiring catflags = 0, to avoid the effect of bad weather conditions (e.g., clouds, moon contamination, and large seeing).We assumed that the time lag between the B-band and ZTF g-band continuum fluxes was much smaller than that of the Hβ emission line, as is the case according to previous continuum reverberation studies (e.g., Jha et al. 2022;Netzer 2022).For example, Wang et al. (2023) reported that a typical size of the continuum emitting region is a factor of ∼8 smaller than that of Hβ.By directly testing the lag between the B and g bands of our sample, we found that the lags between the two broadband light curves were much smaller than the Hβ time lags for the SAMP AGNs (A.K. Mandal et al. 2024, in preparation).By combining the SAMP and ZTF light curves, we obtained improved light curves with better temporal coverage with a cadence of less than 3 days, increasing the constraints of the continuum variability and the reliability of the lag measurements.

Spectroscopy
We carried out the spectroscopic monitoring using two telescopes, the Shane 3 m telescope at the Lick Observatory and the 2.4 m telescope at the MDM Observatory.For the Lick observations, we utilized the Kast double spectrograph,21 which consists of two spectrographs, optimized for red and blue wavelength ranges, respectively.In this study, we focused on the Hβ emission line and used only the red side of the observed spectra.We used the 600 line mm −1 grating, covering the spectral range of 4300-7100 Å, and a dispersion of 2.33 Å pixel −1 until 2016 September.After the upgrade of the CCD in 2016 September, the spectral range was changed to 4450-7280 Å, with a dispersion of 1.27 Å pixel −1 .We used a 4″ slit width to minimize slit loss.Combined with the grating, the Lick spectral setup provided a spectral resolution R of ∼624.We measured the instrumental resolution (FWHM) of 481 km s −1 by utilizing unblended skylines.Calibration frames -i.e., bias, dome flats, and arc lamps (He, Ne, Ar, and Hg-Cd) -were obtained at each night.For most objects, we used a fixed slit position angle (PA) for observations at airmass less than 1.3, while we adopted a parallactic angle for observations at higher airmass.The on-source exposure time of the Lick observations was set between 360 and 1800 s, depending on the magnitude of each target, in order to obtain a signal-tonoise ratio (S/N) per pixel > 15-20 calculated over the entire spectral range.
For the MDM observations, we utilized the volume phase holographic blue grism with a spectral coverage of 3970-6870 Å and a dispersion of 0.715 Å pixel −1 .We used a 3″ slit width before 2017 February, after which we ordered and replaced it with a customized 4″ slit, in order to make a consistent setup compared to the Lick spectroscopy.The corresponding instrumental resolution is R = 617.Calibration frames, including bias, dome flats, and Ar/Xe arc lamps, were obtained at each night.The PA of the slits was set to be the same as the Lick observations.The on-source exposure time of the MDM observations was set between 600 and 2400 s, depending on the magnitudes of the target AGNs.
The SAMP sample covers a large range of R.A., and individual targets show various levels of flux variability.To optimize the monitoring efficiency, we continuously checked the variability and the feasibility of the lag measurements based on the updated light curves on 1-2 month timescales.Note that we can predict the strong variability of the Hβ line emission if we see a strong variability pattern in the photometric continuum light curves in advance.Thus, we reduced the time allocation of relatively nonvarying targets, while we provided more spectroscopic time to promising targets, for which strong variability was detected in the photometric light curves.Consequently, the cadence and the total observed number of epochs varied for each target.In Table 2, we summarize the observations of each object in the final sample.We performed standard spectroscopic data reduction, including overscan subtraction, bias, and flat-fielding using the standard IRAF package.The cosmic-ray rejection was done using the LA-Cosmic task.For the MDM spectra, we used a single sensitivity function, which was averaged over monitoring seasons, because the difference among the various epochs was sufficiently small.In the case of the Lick calibration, we tested the consistency using individual sensitivity functions obtained at each night, which were constructed by fitting the spectra of the spectrophotometric standard stars observed at each night (Oke 1990) with polynomial functions, using a script provided by PypeIt v1.4 (Prochaska et al. 2020a(Prochaska et al. , 2020b)).We found that the results were almost consistent with those from the analysis with IRAF, but provided better consistency at the edge of the spectral range, especially for the 2017-2021 observations.Thus, we adopted the PypeIt reduced spectra for the data obtained over 2017-2021.For the 2016 data, the individual sensitivity function could not be accurately constrained at the blue edge of the spectral range.Thus, we decided to use the reduced spectra based on the IRAF analysis for the 2016 Lick observations.

Flux Recalibration
To reliably measure the spectral variability of the Hβ line emission, we first perform flux recalibration using the nonvarying narrow [O III] line emission.This rescaling approach described by van Groningen & Wanders (1992) is often adopted by various RM studies, due to the nonintrinsic variation caused by the slightly inconsistent spectral resolution and the slit loss due to the nightly change of focus, seeing, and centering of the AGN in the slit, etc. (e.g., Peterson et al. 1995).
We utilize the Python package mapspec22 (Fausnaugh 2017), which follows the same approach proposed by van Groningen & Wanders (1992), but works in a Bayesian framework, enabling an assessment of the uncertainties of the calibration.As We further calibrate the absolute flux level of the spectroscopy by matching the synthetic V-band light curves with the photometric V-band light curves (or R-band for two objects; J1105+671 and PG 2251+113 at z>0.3) using PyCALI.In this process, we obtain an average scale factor and rescale the Hβ light curve for each object.The synthetic V(R)-band flux is calculated by performing synthetic photometry on the mapspec-calibrated spectra.
where N is the number of spectra,  2 to demonstrate the quality of the flux scaling.
We note that for AGNs with a weak [O III] line, this procedure introduces a large systematic error, since it is difficult to define the [O III] line profile because of the blending of strong Fe II and [O III] lines.In our sample, five objectsnamely, J0939+375, PG 1121+422, PG 1322+659, PG 1440 +356, and J1526+275-show a relatively small [O III] equivalent width.As a different approach, we obtain a flux scale factor for each epoch by matching the synthetic V-band magnitude with the V-band magnitude from the interpolated photometry light curve based on the DRW model provided by PyCALI.By comparing the results from two different calibration approaches, we find that only one target, J1526 +275, with the weakest (almost no narrow) [O III] line shows noticeably better cross-correlation results from the photometrybased scaling.Therefore, we adopt the photometry-based scaling result for J1526+275, while we use the [O III]-based scaling results for the rest of the targets.

Mean and rms Spectra
We generate the mean and rms spectra using the fluxrecalibrated spectra of each target, as done by previous studies (e.g., Peterson et al. 2004;Barth et al. 2015; Paper II).The mean spectra are calculated by averaging the flux of all epochs, while the rms spectra are generated using the definition å = -- where N is the number of spectra, and f λ,i and l f are the flux density of the spectrum from the ith epoch and the mean spectrum, respectively.
We provide two sets of rms spectra, using the individual epoch's spectra with and without subtracting continuum and narrow emission lines (see Section 3.3 for the decomposition process).It has been shown that the rms spectra generated without subtracting the continuum and narrow lines can suffer from a bias in the line width measurement, because of the different variations between the continuum and emission-line fluxes (Barth et al. 2015;Wang et al. 2019).As presented in Figure 2, we find a small residual at the location of [O III] λ5007, indicating the high quality of the flux calibration.The noticeable residual of [O III] in the rms spectra of some targets, together with other small features caused by the telluric lines or the CCD cosmetics and cosmic rays, are masked out (as indicated by the orange color in Figure 2) and interpolated before the line width calculation.Note that because we have different wavelength coverages, owing to three instrumental setups of Lick spectroscopy, there are small rises of the continuum at the blue side of Hβ (e.g., PG 1440+356 and Mrk 1501) due to the CCD effect.

Spectral Decomposition and Measurements
We perform spectral decomposition on the mean spectrum to remove the host contamination, if any, and measure the line flux and width of the broad Hβ line.The same procedure is adopted to fit the spectrum of each epoch to derive the light curve of the Hβ line flux.Compared to the linear fit of the continuum under the Hβ line profile, the spectral decomposition approach can better isolate the blended components, e.g., Hβ and Fe II.Thus, the decomposition approach has been adopted in a number of recent studies, particularly those based on high-quality spectra (e.g., Park et al. 2012b;Barth et al. 2015;Hu et al. 2015Hu et al. , 2021;;U et al. 2022).
We adopt the spectral decomposition procedure of the previous studies by Shen et al. (2011Shen et al. ( , 2019) ) with slight modifications.Here, we summarize the basics of the fitting procedure.First, we generated a pseudocontinuum model by combining a power law, an Fe II component, and a host galaxy stellar population model, in order to fit the observed spectra within two line-free windows, i.e., [4450,4600] and [5050,5550] Å in the rest frame.We apply velocity shift and Gaussian velocity broadening to the iron and host templates, respectively.For the Fe II model, we test two commonly used templates provided by Boroson & Green (1992) and Kovačević et al. (2010).In the case of the template of Kovačević et al. (2010), it is difficult to use consistent flux ratios of different emission-line groups for individual epochs, because of the variation of the wavelength coverage due to the different spectroscopy setups of the MDM and Lick observations.Thus, we decide to use the Fe II template of Boroson & Green (1992).However, we find that the Fe II template of Kovačević et al. (2010) provides a more accurate fit for several objects (e.g., Mrk 1501 and Mrk 1014; see Barth et al. 2015;U et al. 2022).Thus, we adopt the fits using the Kovačević et al. (2010) iron template for these objects.In the case of the host galaxy stellar Note.For continuum light curves, the quantity F represents the flux density in units of 10 −15 erg s −1 cm −2 Å −1 , while for Hβ (band = Hβ), the quantity F represents the integrated Hβ flux in units of 10 −15 erg s −1 cm −2 .
(This table is available in its entirety in machine-readable form.) population model, we used a single-burst 10 Gyr old stellar population model with solar metallicity (Bruzual & Charlot 2003).As our targets are moderate-to-high-luminosity AGNs, the host galaxy stellar absorption line features are negligible in the mean spectra of most targets.We present one example of the spectral decomposition using the lowest-luminosity AGN, J1217+333, in our sample in Figure 3.As expected from its low AGN luminosity, the host component, i.e., the strong Mg Ib λλ 5167,5173,5184 triplet absorption line, is clearly detected.In contrast, most of the other objects show no clear signs of host galaxy absorption lines in their mean spectra.
After subtracting the best-fit pseudocontinuum model, we fit the residual emission lines with Gaussian models by accounting for various kinematical features of individual emission lines.Specifically, we use three Gaussians for the broad Hβ line and one Gaussian for the narrow Hβ component.For [O III] λ4959 and [O III] λ5007, we use two Gaussian models for the central and wing components.In addition, one narrow Gaussian model and one broad Gaussian model were used for fitting the narrow and broad components of He II λ4687.All the narrow line centers and widths are tied together.

Hβ Light Curves
We measure the Hβ line flux based on the decomposition and generate the light curves over the 6 yr monitoring period.Note that we apply a window to each epoch's spectrum, within which the Hβ line flux is summed, as listed in Table 3.These additional windows are added to avoid any overfitting to the noise at the edge of the line profile.4): the observed-frame lags of the ICCF centroid (τ cent ), peak (τ peak ), and JAVELIN, respectively, through weighting and the alias removal procedure.The lags and their two uncertainties are determined from the median, while the 16th-50th/84th-50th percentile intervals are determined from the unweighted lag posterior distribution covered by the primary peak.For targets with asterisks, the lags are calculated based on part of the light curves.Columns (5)-( 8): lag reliability parameters: r max represents the maximum correlation coefficient; p r max ( ) is derived from the PyI 2 CCF simulation, which indicates how large the chance is to obtain the observed r max from random light curves; and f peak is the fraction of the posterior distribution within the selected primary peak.Seven objects with less reliable lag measurements indicated by r 0.6 max  , p r 0.2 max  ( ) , or f peak 0.6 are listed at the bottom.
(This table is available in machine-readable form.) Note that while both the Lick and MDM spectra were aligned to the same reference spectrum during the flux recalibration process using the nonvarying [O III] emission line, there can still be a small offset between two sets of light curves from two different telescopes, Lick and MDM, due to slightly different aperture effects (Peterson et al. 1995).We (blue) and JAVELIN τ JAV (orange), as well as the applied weight (dotted-dashed line) for searching the primary peak (see Section 4.1), are presented in the lower right panel.The blue and orange shadowed areas indicate the range of the primary peak for τ cent and τ JAV , respectively.The vertical solid and dashed lines represent the location of the lag and its upper and lower limits, calculated as the median, 16th, and 84th percentiles within the primary peak range, respectively.τ cent is adopted as the final lag measurement, by which the continuum light curve is shifted and matched with the Hβ light curves, as visualized by the gray points in the lower left panel.The two lag reliability indicators, i.e., the r max and p r max ( ) (see Section 4.3 for details), are displayed in the upper right panel.
further match the two light curves from the Lick and MDM telescopes by assuming the spectra obtained at a similar time should have the same flux levels.By searching for pairs of close epochs within one day between the Lick and MDM observations, we calculated the median scaling ratio between the two spectra using these pairs.Finally, the MDM light curves are scaled with these ratios to match with the Lick light curves.The final continuum and Hβ light curves are presented in Table 4.
For each object, we calculate the noise-corrected fractional variability F var of both the continuum and Hβ light curves, as this quantity is frequently adopted to represent the level of variability (e.g., Peterson et al. 2004;Bentz et al. 2009b).The F var and its uncertainty is defined by Rodríguez-Pascual et al. (1997) and

12
The Astrophysical Journal, 962:67 (31pp), 2024 February 10 Woo et al.Edelson et al. (2002) as where N is the number of epochs, 〈f〉 and σ are the mean flux and standard deviation of the light curve, and δ is the rms of the individual flux uncertainties.We list the F var of the continuum and Hβ in Table 2.For our sample, the F var of the continuum ranges from 0.07 to 0.21, while the F var of Hβ ranges from 0.05 to 0.20.Note that 30 out of 32 objects in the sample show continuum F var 0.1, indicating the clear detection of the variability of the majority of the targets.

Hβ Lag Measurement
The emission-line lags have been measured in a number of different ways in previous RM studies.These methods are categorized into two main classes: traditional cross-correlation analysis and the methods with a statistical approach for describing AGN variability, e.g., DRW models (e.g., Kelly et al. 2009;MacLeod et al. 2010).The former class includes the interpolated cross-correlation function (ICCF; e.g., Gaskell & Peterson 1987;Peterson et al. 1998)  comprehensive comparison among the second-class methods is not yet available, these methods generally share a similar basic algorithm and at least some of them provide consistent results (e.g., Grier et al. 2017Grier et al. , 2019)).
We adopt the most commonly used approach, ICCF, as the primary method for the lag measurements of SAMP targets.As the ICCF method has been adopted by the majority of previous studies, it is possible to directly compare with the previously reported measurements.We use the Python package PyCCF (Peterson et al. 1998;Sun et al. 2018) to perform the lag measurements, by calculating the cross-correlation coefficient r between the continuum and Hβ emission-line light curves, after linearly interpolating one light curve to match the time grid of the other light curve.One light curve is shifted by a series of τ values in a searching window and, as a function of τ, the ICCF is calculated and the cross-correlation coefficient r is determined.In this process, either the continuum or the emission-line light curve is interpolated to calculate the ICCF, and the final ICCF is obtained by averaging the two ICCFs.Then the centroid (τ cent ) or the peak (τ peak ) is determined using the range in the averaged ICCF, where r is larger than 80% of its maximum, and adopted as the lag of the two light curves.We consider τ cent as the primary ICCF lag measurement, as done in previous studies (e.g., Peterson et al. 2004), while we also present the τ peak of each AGN for comparison and completeness in Table 5.The lag uncertainties are estimated with the flux randomization/random subset sampling method, which randomizes the flux of each epoch based on its uncertainty and randomly selects a subset of epochs in the light curve for each simulation (Peterson et al. 1998).
In the ICCF analysis, we use a searching window of [−600, 600] days, with a step of 1 day.Note that this window is roughly 30% of our 6 yr baseline (∼2000 days on average) and at least a factor of 3 larger than the expected lag, estimated based on the size-luminosity relation (Bentz et al. 2013).For two targets, PG 1100+772 and PG 2251+113, the searching window is less than a factor of 2 of their expected lags.Thus, we test the effect of the size of the searching window by increasing the searching window as [−1000, 1000] days, finding that the lags are essentially the same.For J1415+483, which is the only target that was monitored for 5 yr, we use a searching window of [−445, 445] days, corresponding to ∼30% of its baseline.In the case of the JAVELIN analysis, we initially use the same searching window.However, we finally use a smaller searching window, [−360, 360] days, because we find JAVELIN tends to provide strong aliases if we use a larger searching window, and sometimes unacceptably large lag values are reported (i.e., >500 days or <−500 days).Notes Column (1): object names.Column (2): the final rest-frame lags.We adopt the ICCF τ cent for all objects except J1526+275, for which we prefer τ JAV as the final lag measurement (see Section 4.3).Six unreliable measurements are listed at the bottom, with one object showing negative lags not displayed.Column (3): extinction-corrected and host-contamination-removed (if any) AGN luminosities in units of 10 44 erg s −1 .Columns (4)-( 5): FWHM and σ line measured from the mean spectrum.Columns (6)-( 7): FWHM and σ line measured from the rms spectrum.Column (8): virial products calculated using the final lags in Column (2) and the rms spectrum σ line in Column (7).Column (8): final BH masses calculated by multiplying the virial products in Column (7) with the virial factor f = 4.47 (Woo et al. 2015).
(This table is available in machine-readable form.) In addition, we test the improvement of the lag measurements by using subsamples of the 6 yr light curves, finding that cross-correlation results are more reliable with higher crosscorrelation coefficients, because of the densely distributed observations in particular seasons.Note that the total 6 yr light curve is composed of six seasons separated by seasonal gaps.Specifically, we find an improvement relative to the total light curves when we use only the 2018 observations for J0939 +375, the 2020-2021 observations for J1203+455, the 2016-2017 observations for PG 1545+210, and the 2019-2021 observations for PG 1427+480, presumably owing to the much higher-quality light curves during these specific seasons.Under such circumstances, we adopt the measurements based on the light curves of these specific seasons.9.The same as Figure 4, but for J1203+455, PG 1202+281, and J1217+333.
Both the posterior distributions of the ICCF and JAVELIN can show multiple peaks, as a feature of sparsely sampled multiple-year data with seasonal gaps (e.g., Grier et al. 2019;Homayouni et al. 2020;Yu et al. 2023), which is likely to be caused by quasiperiodic variations, the mismatch of weakly variable features (Homayouni et al. 2020), and seasonal gaps.Following Grier et al. (2019), we employ an alias identification procedure to remove these aliases and identify the primary peaks for measuring lags.We apply a weight function to the posterior distribution of the ICCF and JAVELIN, which is a convolution of the following two components.The first component is a probability function based on the overlapped fraction, defined as 10.The same as Figure 4, but for VIII Zw218, PG 1322+659, and J1415+483.
where N(τ) and N(0) are the number of overlapped points between the continuum and emission-line light curves with or without shifting one light curve by τ, respectively.This component helps to reduce the weights of the lag values that are similar to seasonal gaps, because no real data are overlapped between the light curves of the continuum and Hβ if one light curve is shifted by these values.The second component is the autocorrelation function (ACF), which represents how fast the continuum variability is.If the continuum light curve varies slowly, then the distribution of the ACF is wide, indicating that the seasonal gaps are less likely to affect the lag measurement.On the other hand, fast continuum variability leads to a narrow ACF distribution, and in this case, the important features of the variability pattern are more likely blocked by the seasonal gaps (Homayouni et al. 2020).Thus, we combine the overlap probability function (i.e., Equation ( 6)) and the ACF to generate a weight function, which then is multiplied by the lag posterior distribution of the ICCF and JAVELIN.
To identify the primary peaks in the weighted posterior distribution, we first smooth the distribution by a Gaussian kernel with a width of 12 days, which was determined based on experiments and visual inspection, as similarly adopted by previous studies (e.g., Grier et al. 2019;Homayouni et al. 2020;Yu et al. 2023).By identifying the highest peak as the primary peak in the smoothed distribution, we define its range between the two adjacent local minima.Then we remove the posteriors outside this range and measure the lag from the truncated distribution.Note that we use the weighted posterior distribution for defining the range of truncation.Finally, we determine the lag and its uncertainty by measuring the median and 16thto-84th quadrature of the truncated unweighted posterior distribution.
Based on the aforementioned analysis, we measure the Hβ lag of the sample (see

Effect of Detrending
A long-term trend in the line light curves can bias the crosscorrelation result because the long-timescale (low-frequency) trend can overwhelm the relatively shorter intrinsic pattern due to the lag (e.g., Welsh 1999;Denney et al. 2010;Zhang et al. 2019).In such cases, the detrending of the light curves may improve the lag analysis.
We examine the effects of long-term trends by detrending the continuum and/or Hβ light curves using a first-order polynomial, which is usually enough to remove the various linear trends between the continuum and emission-line light curves.Higher-order polynomials are more dangerous to use, since they are suspected to introduce artificial signals into the light curves and largely reduce the cross-correlation strength at the same time (Peterson et al. 1995).We find that only for two objects, PG 1100+772 and PG 1202+281, r max is significantly improved after detrending.Thus, we decide to use the results based on detrending for the two objects.For completeness, we present the lag analysis for these two objects without detrending in Appendix A.
Based on the optical and radio variability of a radio-loud AGN, 3C 273, Li et al. (2020) suggested that the underlying physics of the detrending process is to correct for the contamination of jet emission in the optical continuum light curves, since jet emission has no corresponding effect on the emission-line light curves.We note that PG 1100+772 is a radio-loud AGN with a radio loudness R = f 5GHz /f 2500 ∼ 400, suggesting that PG 1100+772 can be treated as an analog of 3C 273 (Zhang et al. 2019).However, it is difficult to disentangle the jet contribution from the observed continuum light curves for this target, due to the lack of radio monitoring data during our reverberation campaign.It is possible that radio-loud AGNs may introduce scatter in the size-luminosity relation if this scenario is correct.A larger sample of radio-loud AGNs is required in RM studies to verify this scenario.

Lag Reliability
In this section, we investigate the reliability of the lag measurements from the cross-correlation analysis.The two main criteria of lag reliability are the goodness of the crosscorrelation between the continuum and emission-line light curves, and the constraints against artificial signatures in the light curves, i.e., aliases.In general, the lag measurements are more reliable if the input light curves have a higher cadence, more regular sampling, a longer time baseline, and stronger nonquasiperiodic variability patterns.However, well-defined quantitative criteria are not yet available, while a few new attempts have been reported in the literature.
One of the criteria of reliability is the maximum crosscorrelation coefficient r max , which is the maximum value in the ICCF, representing the strength of the correlation between two light curves.Previous studies have often used a threshold of r max > 0.6 or 0.5, to adopt reliable lag measurements (e.g., Grier et al. 2017;U et al. 2022).However, r max alone is not sufficient to determine the lag reliability.For instance, r max would be relatively high even though the light curves have very few data points, but show strong linear variation, or even if the variability amplitude is small compared to the flux uncertainties, but some features in the light curves are well matched.In these cases, the measured lag can be biased or reflect randomly uncorrelated light curves.Not to mention that flux uncertainties are not taken into account in the r max calculation.
Recently, a new method PyI 2 CCF23 has been proposed to assess the lag reliability (Guo & Barth 2021;U et al. 2022;H. Guo et al. 2024, in preparation), using a null hypothesis to evaluate the probability that the observed cross-correlation can be equally obtained by two uncorrelated red-noise light curves.This idea originates from X-ray reverberation studies and has been used to assess the correlation of multiwavelength variability of AGNs (Uttley et al. 2003;Arévalo et al. 2008;Chatterjee et al. 2008).Using the publicly available PyI 2 CCF, we investigate the reliability of our lag measurements.We generate 10 3 realizations of a pair of mock light curves for each object by keeping the same cadence and sampling as the observed light curves using DRW models.Then, by counting the number of positive lags (τ > 0) with the r max value higher than the observed r max among all simulations, we determine the probability t> p r max 0 ( ) as a lag significance indicator.In this process, we keep the aforementioned searching window as defined in Section 4.1.
We further examine the fraction of the primary peak ( f peak ) in the posterior distribution that describes how much of the posteriors are within the range of the primary peak.For example, an f peak of less than 0.6 is often considered as indicating that there are possible significant solutions other than the identified primary peak (Grier et al. 2019;Homayouni et al. 2020;Yu et al. 2023;Malik et al. 2023).We list the f peak of τ cent of all targets in Table 5.
We present the three criteria of lag reliability-r max based on observed light curves, p r max ( ) based on simulations, and f peak calculated from the posterior distribution-for individual AGNs in Table 5.Note that for the majority of our targets, the measured lag is reliable based on these three assessments.
In Figure 15, we directly compare the two indicators, i.e., r max and p r max ( ), along with the color-coded ( f peak ).We quantitatively define the best lag measurements, by requiring r 0.6 max  , p r 0.2 max  ( ) , and f peak 0.6.Note that there is no strict cutoff of p r max ( ) between the best and less reliable lag measurements.We adopt t> p r 0.2 max 0  ( ) following the previous study by U et al. (2022).In these assessments, we find that the measured lags of seven targets-namely, J0801 +512, J1059+665, J1105+671, J1203+455, PG 1545+210, J1935+531, and PG 2251+113-do not satisfy the lag reliability criteria.By excluding these seven targets, we obtain the best Hβ lag measurements for 25 AGNs.In the case of PG 0026+129, the obtained lag is consistent with zero within 1σ uncertainty, suggesting that the lag is not resolved.
Finally, we perform a visual inspection of all light curves to check how well the continuum light curve matches the emission-line light curve after shifting the continuum light curve by the measured lag.This method is qualitative, providing an additional check on the reliability of the measured lag.We present the two light curves, after shifting by the lag and scaling with the median and standard deviation of the fluxes, in the lower panel of each target in Figure 4. Overall, we find a good match between the shifted continuum and the Hβ light curves of the majority of targets in the sample, confirming the reliability of the lag measurements.
As a consistency check, we compare the ICCF and JAVELIN measurements in Figure 16, finding that they are generally consistent.However, six objects-namely, Mrk 1501, J1026+523, J1059+665, J1105+671, PG 1322 +659, and PG 2251+113-show a large discrepancy, with the JAVELIN lag larger than the ICCF lag by 150 days or smaller by −150 days while the ICCF lags are within [0, 100] days.Three of them are already identified as less reliable, while for the other three AGNs the JAVELIN lag seems overestimated or underestimated, because of the combined effect of seasonal gaps and quasiperiodic variability, based on the visual inspection of their light curves.This result suggests that although the typical error of the JAVELIN lag is indeed smaller than the ICCF one, the JAVELIN measurements are sometimes more easily affected by seasonal gaps than the ICCF ones and generate spurious measurements, particularly for sparsely sampled light curves.Thus, we decided not to disqualify these ICCF lags as unreliable, despite the discrepancy with the JAVELIN results.
A simple statistical approach to quantifying the lag reliability of a given sample is to examine the sample's false-positive rate, i.e., what fraction of lag measurements are false detections caused by inopportune features, seasonal gaps, as well as inappropriately identified primary peaks and alias identifications.The incidence of negative lags over the entire sample can be an indicator of the false-positive rate (e.g., Shen et al. 2016;Grier et al. 2017Grier et al. , 2019)).In our SAMP sample, there is one object with a negative τ cent and one object with a negative τ JAV .
If we conservatively count these two measurements, the falsepositive rate is ∼6% of the SAMP sample, indicating the high quality of our lag measurements.

BH Mass and Eddington Ratio
In this section, we present the M BH determination of the sample, by measuring the width of the Hβ emission line.We also determine the Eddington ratio of individual AGNs using continuum luminosity.We compare the distribution of the M BH , luminosity, and Eddington ratio of the SAMP AGNs with that of the previous reverberation-mapped AGNs.

Hβ Line Width Measurements
We measure the Hβ line width from the decomposed broad Hβ emission-line profile using the mean spectra as well as the broadline-only rms spectra.As a velocity measure, we adopt the FWHM and the line dispersion σ line (the second moment of the line profile), following the definition by Peterson et al. (2004).Note that we conservatively define two windows to represent the local continuum in the rms spectra for line width measurements, as shown in Figure 2. Using the measured instrumental resolution of 481 km s −1 from skylines, we correct for the instrumental broadening by subtracting the instrumental resolution from the measured widths in quadrature.Since the spectral resolution is measured from skylines, which uniformly fill the 4″ slit, the actual resolution for a point source is likely to be smaller.However, this effect is insignificant, because Hβ lines are intrinsically very broad.Since the change of the instrumental broadening during our observations is relatively small compared to the broadline width, we adopt a single representative spectral resolution, 481 km s −1 .

M BH Measurements
We determine the M BH based on Equation (1), using the τ cent from the ICCF analysis as the lag and the line dispersion of Hβ measured from the rms spectra (σ rms ) as the velocity, along with the average f = 4.47 ± 0.43 (Woo et al. 2015).Note that while we adopt a single f value in this work, we will present the velocity-resolved lag measurements and dynamical modeling to derive the f factor for each AGN in the future.
For determining M BH , the line dispersion (σ rms ) is commonly used in RM studies, because it best recovers the virial relation (Peterson et al. 2004;Park et al. 2012b).It is also suggested that the line dispersion is less sensitive to other factors, e.g., orientation, compared to the FWHM (Collin et al. 2006;Wang et al. 2019).On the other hand, the FWHM is frequently used in single-epoch M BH estimation, as it is easier to measure than the line dispersion, based on low-S/N spectra.We provide the M BH of the sample based on the line dispersion in Table 6, but the M BH with FWHM can be easily derived using Table 6.
We also estimate the bolometric luminosity, using the monochromatic luminosity L 5100 along with a constant bolometric correction factor κ = 9.26 (e.g., Richards et al. 2006;Shen et al. 2011).For calculating the Eddington luminosity L Edd , we assume L Edd = 1.3 × 10 38 M BH /M e .
In Figure 17, we present the M BH and bolometric luminosity of the sample, in comparison with those of the previously reverberation-mapped AGNs, by combining the collection of Dalla Bontà et al. (2020) and the AGNs from SEAMBHs (Du & Wang 2019;Hu et al. 2021;Li et al. 2021).The M BH of the SAMP AGNs range from 10 7 to 10 9 M e and the Eddington ratio ranges from 0.05 to 1.82, with a median of 0.17, which is higher than that of the sample of Dalla Bontà et al. (2020), but lower than the AGNs in SEAMBHs.

Hβ BLR Size-Luminosity Relation
As the number of Hβ lag measurements has increased over the last decade, it has been demonstrated that the Hβ sizeluminosity relation has a considerably larger scatter and a shallower slope than was previously accepted, indicating a more complex nature of the relation (e.g., Du et al. 2016Du et al. , 2018a;;Grier et al. 2017;Du & Wang 2019;Martínez-Aldama et al. 2019;Hu et al. 2021;Li et al. 2021).In particular, it is claimed that super-Eddington AGNs tend to be systematically offset from the relation of lower-Eddington-ratio AGNs (Du et al. 2016;Li et al. 2021).In this section, we investigate the Hβ BLR size-luminosity relation by combining our new measurements of high-luminosity AGNs with previous measurements from the literature.

Best-fit Slope and Scatter
We utilize the collection by Dalla Bontà et al. (2020), who compiled a sample of AGNs with Hβ lag and L 5100 after correcting for host galaxy contamination.As this collection does not include a large fraction of high-luminosity AGNs and AGNs without available HST images, we collect a number of AGNs from various other sources, i.e., the SEAMBHs project (e.g., Du et al. 2016Du et al. , 2018a;;Hu et al. 2021;Li et al. 2021), the SDSS-RM project (Grier et al. 2017), LAMP 2016 (U et al 2022), and other recent studies (Hu et al. 2021;Li et al. 2021;Malik et al. 2023).The total sample consists of 242 Hβ Figure 15.Assessment of the lag reliability using the three criteria, i.e., the maximum cross-correlation coefficient (r max ), the probability that uncorrelated light curves would produce a cross-correlation coefficient of at least r max (p r max ( )), and the fraction of the primary peak in the posterior distribution ( f peak ).The best lags are defined with r max 0.6 (vertical dashed line), p r max ( ) < 0.2 (horizontal dashed line), and f peak > 0.6.and p r max ( ) < 0.2, but with f peak 0.6.lag measurements (including multiple measurements of NGC 5548), among which 30 AGNs are based on the SAMP results.We exclude two SAMP AGNs, namely PG 0026+129 and J1935+531, for which the obtained lags were negative or not larger than zero by 1σ uncertainty (see Table 5).Note that by adding these two AGNs, we find an insignificant change of the slope and scatter.
We perform a linear regression to derive the best-fit relation as where the Hβ BLR size is in light-days and the monochromatic luminosity λL 5100 is expressed in units of 10 44 erg s −1 .The α and K are the slope and intercept.To account for the asymmetric uncertainties of the lags, we performed MCMC regression fits using the Python package emcee24 and adjusted the likelihood function, so that it uses the upper error when the model value is larger than the data and uses the lower error in the opposite case.In brief, the regression minimizes the following quantity: where x err,i and y err,i are the errors of variable x i and y i , respectively, x 0 is the normalization x, and σ int is the intrinsic scatter.
We obtain a best-fit slope of -+ 0.402 0.022 0.020 using the total sample of 240 AGNs with an intrinsic scatter of 0.232 dex (Figure 18, left).As we find more similarity between SAMP and LAMP in terms of spectroscopic observations, spectral analysis, and lag measurements, we use a combined sample of the 24 best SAMP measurements (after excluding PG 0026 +129 due to its unresolved lag) and 23 LAMP measurements, finding a best-fit slope of -+ 0.444 0.035 0.036 with an intrinsic scatter of 0.177 dex (Figure 18, right, and Table 7).Note that most of the SAMP objects lie below the previous relation (i.e., α = 0.533), which was defined based on ∼40 AGNs (∼70 lag measurements) by Bentz et al. (2013).
It is unlikely that the deviation of the SAMP AGNs, with the observed-frame lag up to 254 days, is due to the underestimation of lags, since the SAMP results are based on the 6 yr (∼2000 days) monitoring data with a ± 600 days lag search window.SAMP AGNs are generally luminous (i.e., L 5100 > 10 44 erg s −1 ) and the host galaxy contamination is not significant, as demonstrated by the lack of stellar absorption lines in the spectra.Based on the spectral decomposition, we measure the host fraction at 5100 Å.The average host fraction of our sample is 0.07, ranging from 0 to 0.28.Thus, the continuum luminosity is unlikely to be overestimated because of the host contribution.
Our results suggest that the slope is likely to be shallower than the popularly used 0.533 slope of Bentz et al. (2013) and that if the previous relation is utilized, single-epoch M BH , particularly for high-luminosity AGNs at high z, would be overestimated.
Note that these Hβ lag measurements collected from the literature were not consistently determined, since various studies adopted a number of different methods in their analysis.Hu et al. 2021;Li et al. 2021; brown triangles).The three dashed lines indicate the Eddington ratios from 0.01, to 0.1, to 1.0, as labeled to the left above the lines.Our SAMP final sample has moderately high Eddington ratios.Bottom: the distribution of the Eddington ratios for SAMP and other samples.Thus, there could be various systematic uncertainties in deriving the size-luminosity relation based on this heterogeneous sample.Nevertheless, we note that the best-fit relation based on the total sample is relatively tight, with an rms scatter of ∼0.3 dex and an intrinsic scatter of 0.234 dex, without requiring two different relations, respectively, for sub-Eddington and super-Eddington AGNs (see Li et al. 2021).It is important to perform a consistent and uniform analysis of the cross-correlation and error measurement, in order to derive a better size-luminosity relation.We will revisit the sizeluminosity relation with uniformly measured lags and investigate the systematic effects of AGN parameters on the sizeluminosity relation in our next study.

Deviation from the Size-Luminosity Relation
We investigate whether the deviation from the sizeluminosity relation correlates with any key parameters of AGNs.Using the best-fit relation obtained based on all targets (Case 1), we calculate the deviation of each target and compare with the Eddington ratio in Figure 19 (right).It was previously reported that super-Eddington AGNs tend to have smaller Hβ sizes for a given AGN luminosity and that this deviation is stronger with higher Eddington ratios or high mass accretion rates.For example, Du et al. (2015) claimed a negative correlation between the mass accretion rate and the deviation.However, this could be naturally caused by the zone of avoidance, since there are typically no type 1 AGNs with the  Note that the decreasing trend with the Eddington ratio is partly due to the zone of avoidance (no type 1 AGNs with FWHM <1000 km s −1 ) and self-correlation.The dotted line indicates the upper envelope of the distribution defined by AGNs with FWHM=1000 km s −1 for a fixed luminosity, L 5100 = 10 44.5 and 10 45 , respectively.
FWHM of the broad Hβ line being less than 1000 km s −1 (see the dashed line in Figure 19).For example, for a fixed AGN luminosity (i.e., L 5100 = 10 45 erg s −1 ) in the Eddington ratio between ∼0.3 and 1, there is an upper envelope above which no AGN can be located, because the widths of the Hβ broad line are not typically smaller than 1000 km s −1 .Note that this effect is not strong at lower Eddington ratios (i.e., < ∼0.3), leading to a negative correlation only in the high-Eddington regime.
The decreasing trend, particularly at the high Eddington ratio, is naturally expected (Figure 19, right) due to the selfcorrelation between the Eddington ratios and the deviation, since the former scales with L/R and the latter scales with R/L 0.5 , as pointed out by Fonseca Alvarez et al. (2020).In other words, if the measured R BLR is somewhat smaller than expected by the best fit, then the M BH would be smaller.Consequently, the Eddington ratio of these targets would be systematically larger.Thus, the decreasing trend between the deviation and the Eddington ratio is expected, due to the selfcorrelation.
To overcome these artificial trends due to the zone of avoidance and the self-correlation, we use the strength of the Fe II emission (R Fe ), which is the ratio between the Fe II emission flux in the range of 4434-4684 Å and the Hβ emission flux, for comparing with the deviation (Figure 19, right).As R Fe is closely related to Eigenvector 1 in principal component analyses of AGN properties, which is typically interpreted as an accretion rate parameter, we use R Fe as an indicator of the Eddington ratio.Note that for this practice, we only use the subsample including SAMP, LAMP, and other AGNs with available R Fe in the literature.However, we find no strong trend between the deviation from the size-luminosity relation and R Fe , suggesting no systematic effect of the accretion rate on the deviation from the BLR size.In contrast, the SEAMBHs project reported that the deviation of their sample correlated with a dimensionless accretion rate (Du et al. 2014(Du et al. , 2015) ) as well as the FWHM/σ ratio of the Hβ line and R Fe (Du et al. 2016;Du & Wang 2019).Note that we reproduce the same result as Du et al. (2016) and Du & Wang (2019) using their sample.However, by adding more luminous sub-Eddington AGNs (e.g., SAMP AGNs), their correlation becomes weaker.We also note that by adopting a shallower slope of 0.4 in the size-luminosity relation, we obtain a systematically smaller correlation coefficient between the deviation and either R Fe or Eddington ratio.A full analysis based on a large sample of reverberation-mapped AGNs with consistently measured Hβ lag, line width, M BH , and R Fe is required to clearly investigate the dependency of the sizeluminosity relation on the Eddington ratio.We will present more detailed results based on our reanalysis of the Hβ lag measurements and the size-luminosity relation using a large archival sample (S.Wang et al. 2024, in preparation).In this section, we discuss several scenarios to explain the observed slope of the Hβ BLR size-optical luminosity relation.First, we expect that the BLR size is proportional to the 0.5 power of the ionizing luminosity if the ionization parameter (U) and hydrogen gas density (n H ) are similar for all AGNs, since U is proportional to L ion /(4πR 2 n H ), where L ion is the photoionizing luminosity.However, this assumption may not be true, as the gas clouds in the BLR could have a range of U and n H (e.g., Baldwin et al. 1995).If higher-luminosity AGNs have an average higher value of the product, U × n H , then the size of the BLR gets smaller than expected, leading to a shallower slope.Currently, we do not have clear observational evidence for this trend.
Second, the optical luminosity measured at 5100 Å may not properly represent the ionizing luminosity, although the UV continuum luminosity at around the Lyman edge may show a better correlation with a 0.5 slope (see the discussions by Czerny et al. 2019;Fonseca Alvarez et al. 2020).In this scenario, the UV-to-optical flux ratio has to be systematically lower for higher-luminosity AGNs, in order to be consistent with the observation that L 5100 is larger than expected from the 0.5 slope for a given BLR size.Some of the thin-disk models showed that the UV-to-optical flux ratio decreases with increasing bolometric luminosity, decreasing Eddington ratio, and decreasing black hole spin (e.g., Davis & Laor 2011;Castelló-Mor et al. 2016;Czerny et al. 2019).However, model predictions show various trends of the UV-to-optical flux ratio, depending on the model assumptions, i.e., radiative efficiency and wind (e.g., Laor & Davis 2014).In general, the UV-optical spectral energy distribution (SED) in the thin-disk models is somewhat inconsistent with observations (e.g., Davis et al. 2007;Davis & Laor 2011).Clearly, more detailed studies are required to understand the effect on the spectral slope.
As an empirical test, we compare the available Hβ line luminosity with L 5100 for a subsample of our collected Hβ reverberation-mapped AGNs, finding a slightly sublinear relation, L Hβ ∝(L 5100 ) 0.92±0.01(Figure 20).If we assume L Hβ is proportional to the ionizing luminosity, we expect a ∼0.46 slope in the BLR size-optical luminosity relation.These results suggest that the systematic change of the SED slope between the UV and optical ranges may contribute to the deviation from the 0.5 slope by increasing L 5100 for a given BLR size and ionizing luminosity.
Third, super-Eddington AGNs may suffer a shortening of the BLR size, owing to the self-shadowing effect of the slim disk, as detailed by Wang et al. (2014).Gas clouds in the BLR with a high inclination angle (with respect to the rotation axis) receive less ionizing photons than gas clouds with a lower inclination angle, due to the shadowing of the funnel in the inner disk, and this effect in super-Eddington AGNs leads to a shortening of the BLR size.Thus, super-Eddington AGNs have systematically smaller BLR sizes compared to sub-Eddington AGNs.Fourth, the BLR size-luminosity relation is driven by the outer boundary of the BLR, which is defined by the dust sublimation radius at the inner edge of the torus (Suganuma et al. 2006).Interestingly, the dust size-optical luminosity relation also shows a shallower slope than 0.5.For example, Minezaki et al. (2019) reported a 0.424 slope between the K-band torus size and V-band continuum, and Mandal et al. (2024) found a ∼0.4 slope between the torus size based on the WISE W1 band and optical luminosity.Perhaps this is also due to the selfshadowing effect of the slim-disk model, as pointed out by Chen et al. (2023).
We consider a sample selection effect of our collected AGNs.While we try to increase the number of high-luminosity AGNs in this study, there is also a lack of high-Eddington-ratio AGNs in the low-luminosity range.To demonstrate the difference in the Eddington ratio distribution as a function of luminosity, we divide the sample into four luminosity bins and present the Eddington ratio distribution in Figure 21.We clearly notice that the median Eddington ratio is increasing for higher luminosity bins.This trend can naturally cause the shallower slope if higher-Eddington-ratio AGNs tend to have smaller BLR sizes for a given luminosity.In other words, a higher fraction of high-Eddington-ratio AGNs in higherluminosity bins can generate a shallower slope.
We have discussed various effects that could be responsible for the shallower slope of the correlation between the BLR size and optical luminosity.It is puzzling to observe a shortening of the BLR size for a given object.For example, the best-studied low-Eddington AGN, NGC 5548, showed a factor of 5-10 smaller BLR size for a given (or similar) optical luminosity (see Figure 13 in Pei et al. 2017).This suggests that the BLR size can change significantly, without changing the luminosity or Eddington ratio.More detailed studies are required to better understand the scatter and slope of the BLR size-luminosity relation.

Summary
We present Hβ RM results based on the 6 yr (2015-2021) data from SAMP.For a sample of 32 high-luminosity AGNs (L 5100,AGN > 10 44.1-45.6 erg s −1 ) at z <∼0.4,we measure the lags between the Hβ and continuum light curves, using both the ICCF and JAVELIN methods.By applying three reliability parameters (r max , p r max ( ), and f peak ), we quantitatively access the lag measurements and use the accepted measurements to investigate the BLR size-luminosity relation.Our main conclusions are: 1.Among the lag measurements of 32 targets, we report the 25 best Hβ lag measurements, which satisfy r max > 0.6, p r max ( )< 0.2, and f peak <0.6.These new measurements significantly increase the current RM sample at the highluminosity end.2. We compare the JAVELIN lag τ JAV with the ICCF lag τ cent , finding that they are generally consistent with each other, but τ JAV tends to systematically offset due to more aliases.Thus, we adopt the ICCF τ cent results as our final lag measurements.3.By comparing the Hβ BLR size and AGN continuum luminosity, we find that most of the SAMP AGNs are located below the previous relation measured by Bentz et al. (2013), suggesting that the slope is shallower than that expected from a simple photoionization model.We find the best slope of 0.39-0.46by combining with previous Hβ lag measurements in the literature.This result indicates that the single-epoch M BH estimates based on the previous sizeluminosity relation can be overestimated.4. It is possible that the deviation from the size-luminosity relation correlates with the Eddington ratio.However, we do not clearly confirm the correlation, except for that caused by the self-correlation between the AGN luminosity and Eddington ratio.Nevertheless, we detect a hint of correlation, using the Fe II relative strength.A consistent analysis of Hβ lag and error measurements based on a uniform method for the large RM sample is necessary to unveil the nature of the size-luminosity relation.
While the sample size of the reverberation-mapped AGNs has significantly increased over the last decade, due to various studies, including SAMP, there is still a scarcity of very-highluminosity AGNs at L 5100 ∼ 10 46 , and a uniform analysis including quantitative assessment of the measured lag is required to properly investigate the Hβ BLR size-luminosity relation.These issues are beyond the scope of this paper and we will revisit them in the future.

Figure 2 .
Figure 2. The Hβ emission-line profiles in the mean (upper panel) and rms spectra (lower panel) of each AGN.The red and gray lines represent the rms spectra, calculated from the individual spectra with and without subtracting the continuum and narrow lines, respectively.The two vertical dashed lines indicate the peaks of Hβ and [O III] λ5007.Some features caused by skylines or [O III] residuals are masked (orange dashed lines).

Figure 2 .
Figure 2. (Continued.) mapspec subtracts a linear continuum to extract the [O III] line flux within a user-defined window, we first transfer each spectrum into the rest frame in order to use essentially the same [O III] extraction window for all objects.The extraction window of [O III] is defined as [4968, 5055] Å and the two adjacent continuum windows are selected to be [4963, 4968] and [5055, 5060] Å.After extracting the [O III] of each epoch, mapspec matches the [O III] line profile of each epoch with that of the reference epoch using three parameters: a wavelength shift factor, a flux-scaling factor, and a linebroadening factor (i.e., the width of a Gauss-Hermite broadening kernel; Fausnaugh 2017).We test various choices of the [O III] window and the broadening kernel, finding that the results are essentially the same.As a reference epoch, we choose the broadest [O III] profile with S/N > 20 and degrade the [O III] line profiles of the other epochs to align them with the reference.The selected reference epoch is typically one of the epochs with bad seeing conditions and suffers from slit loss of the flux.One of the advantages of mapspec is that the uncertainty of each parameter can be estimated using the half-amplitude of the 16th-84th quadrature interval of the distribution of Markov Chain Monte Carlo (MCMC) samples.For each individual epoch, the derived uncertainty of the multiplicative factor is added in quadrature to the Hβ flux uncertainty.Thus, the epochs with more uncertain [O III]-based flux calibration have larger uncertainties in the Hβ light curves.

Figure 3 .
Figure 3. Decomposition of the mean spectrum for J1217+333 as an example.The black and red solid lines represent the mean spectrum and best-fit model.The continuum model consists of a power-law component (magenta dashed line), an Fe II template (green dashed line), and a host galaxy template (brown dashed line).Broad and narrow Gaussians are displayed using blue and cyan dashed lines, while the total broad Hβ and He II profiles are shown by the orange and purple solid lines, respectively.
F i and δ i are the [O III] flux and flux uncertainty of epoch i, and 〈F〉 is the average of the [O III] flux.The normalized excess variance is the fractional residual scatter of the [O III] line flux, representing the systematic uncertainty of the flux calibration.While several previous studies uniformly added a single value of s nx 2 to each epoch's flux uncertainty in quadrature (e.g., Barth et al. 2015; U et al 2022), we add the uncertainty of the [O III] scaling factor derived by mapspec to the [O III] flux uncertainty in quadrature for each epoch.For completeness, we present the σ nx of each target in Table

Figure 4 .
Figure 4. Light curves of the continuum (upper left) and Hβ (lower left) in units of 10 −15 erg s −1 cm −2 Å −1 and 10 −15 erg s −1 cm −2 , respectively, along with the ICCF (upper right) and the posterior distribution in the observed frame (lower right) for each AGN.The unweighted posterior probability distribution of the ICCF τ cent (blue) and JAVELIN τ JAV (orange), as well as the applied weight (dotted-dashed line) for searching the primary peak (see Section 4.1), are presented in the lower right panel.The blue and orange shadowed areas indicate the range of the primary peak for τ cent and τ JAV , respectively.The vertical solid and dashed lines represent the location of the lag and its upper and lower limits, calculated as the median, 16th, and 84th percentiles within the primary peak range, respectively.τ cent is adopted as the final lag measurement, by which the continuum light curve is shifted and matched with the Hβ light curves, as visualized by the gray points in the lower left panel.The two lag reliability indicators, i.e., the r max and p r max ( ) (see Section 4.3 for details), are displayed in the upper right panel.
, the discrete correlation function (DCF; Edelson & Krolik 1988), and the z-transformed discrete correlation function (zDCF; Alexander 2013), while the second class includes JAVELIN (Zu et al. 2011), CREAM (Starkey et al. 2016), and MICA (Li et al. 2016).It has been demonstrated that for sparsely sampling light curves, DCF and zDCF are less efficient in recovering the lag than ICCF and JAVELIN (White & Peterson 1994; Li & Shen 2019).While a

Figure 16 .
Figure 16.Comparison of the lag measurements in the observed frame between τ cent (ICCF) and τ JAV (JAVELIN).The best and less reliable lag measurements are denoted with the filled and open circles, respectively, using the threshold of = r 0.6 max

Figure 17 .
Figure 17.Top: the distribution of the SAMP final sample (red circle) on the BH mass-bolometric luminosity plane, overplotted with the Dalla Bontà et al. (2020) collection (black dots) and the SEAMBHs AGNs (Du & Wang 2019;Hu et al. 2021;Li et al. 2021; brown triangles).The three dashed lines indicate the Eddington ratios from 0.01, to 0.1, to 1.0, as labeled to the left above the lines.Our SAMP final sample has moderately high Eddington ratios.Bottom: the distribution of the Eddington ratios for SAMP and other samples.

Figure 18 .
Figure 18.Left: Hβ BLR size-luminosity relation of the combined sample of SAMP (red filled and open circles) and literature measurements, including LAMP 2008 and 2016 (Bentz et al. 2009b; U et al. 2022; blue circles), 12 PG quasars from Kaspi et al. (2000) with the optical luminosity from Bentz et al. (2013; blue open squares), SEAMBHs (Du et al. 2016, 2018a; Zhang et al. 2019; Hu et al. 2021; Li et al. 2021; brown open square diamonds), OzDES (Malik et al. 2023; purple open squares), SDSS-RM (Grier et al. 2017; green open circles), and others in the collection of Dalla Bontà et al. (2020; black open circles).The brown dashed line and the light orange lines represent the best-fit relation and 50 realizations randomly drawn from the MCMC chains.As a comparison, we denote the best-fit slope of 0.533 from Bentz et al. (2013; gray dotted-dashed line).Right: Hβ BLR size-luminosity relation based on the 24 best measurements from SAMP (red circles) and 23 measurements from LAMP (blue circles).

Figure 19 .
Figure19.Hβ comparison of the deviation from the size-luminosity relation with the Eddington ratio (left) and strength of Fe II (right).Note that the decreasing trend with the Eddington ratio is partly due to the zone of avoidance (no type 1 AGNs with FWHM <1000 km s −1 ) and self-correlation.The dotted line indicates the upper envelope of the distribution defined by AGNs with FWHM=1000 km s −1 for a fixed luminosity, L 5100 = 10 44.5 and 10 45 , respectively.

Figure 20 .
Figure 20.Comparison of the optical luminosity at 5100 Å with the Hβ line luminosity for a subsample of the Hβ reverberation-mapped AGNs.

Figure 21 .
Figure21.Eddington ratio distribution in each luminosity bin.It is clearly shown that higher-luminosity bins contain on average higher-Eddingtonratio AGNs.

Table 1
Sample Properties and Observation Parameters

Table 3
Hβ Line Flux Extraction Window and Average Flux To verify the quality of the flux calibration, we calculate the normalized excess variance s nx 2 of the [O III] line flux (Barth et al. 2015) using the calibrated spectra of each target: Notes.Column (1): object name.Column (2): extraction window for the line flux measurements in the rest frame in units of Å. Column (3): average flux ( b F H ) of the Hβ light curve as well as its standard deviation in units of 10 −15 erg s −1 cm −2 .(Thistable is available in machine-readable form.)

Table 4
Continuum and Hβ Light Curves

Table 6
Final Rest-frame Lags, AGN Luminosities, Line Widths, Virial Products, and Black Hole Masses Table 5 and Figures 4-14).As a consistency check, we compare the new measurements with the measurements presented in Paper II, which reported the initial Hβ lag measurements of two AGNs, J0801+512 (2MASS J10261389+5237510) and J1619+501, based on the first 3 yr data.The lags measured in Paper II are