The Lick AGN Monitoring Project 2016: Velocity-resolved Hβ Lags in Luminous Seyfert Galaxies

We carried out spectroscopic monitoring of 21 low-redshift Seyfert 1 galaxies using the Kast double spectrograph on the 3 m Shane telescope at Lick Observatory from 2016 April to 2017 May. Targeting active galactic nuclei (AGNs) with luminosities of λ L λ (5100 Å) ≈ 1044 erg s−1 and predicted Hβ lags of ∼20–30 days or black hole masses of 107–108.5 M ⊙, our campaign probes luminosity-dependent trends in broad-line region (BLR) structure and dynamics as well as to improve calibrations for single-epoch estimates of quasar black hole masses. Here we present the first results from the campaign, including Hβ emission-line light curves, integrated Hβ lag times (8–30 days) measured against V-band continuum light curves, velocity-resolved reverberation lags, line widths of the broad Hβ components, and virial black hole mass estimates (107.1–108.1 M ⊙). Our results add significantly to the number of existing velocity-resolved lag measurements and reveal a diversity of BLR gas kinematics at moderately high AGN luminosities. AGN continuum luminosity appears not to be correlated with the type of kinematics that its BLR gas may exhibit. Follow-up direct modeling of this data set will elucidate the detailed kinematics and provide robust dynamical black hole masses for several objects in this sample.


INTRODUCTION
It has been known for about two decades that black hole (BH) mass exhibits a tight correlation with the stellar velocity dispersion of the galactic bulge within which it resides, and this relation holds over several orders of magnitude in black hole mass (e.g., Ferrarese & Merritt 2000;Gebhardt et al. 2000).This relation suggests that supermassive black holes and their host galaxies are tightly linked throughout their lifetimes.A fundamental understanding of the growth of supermassive black holes and their interaction with their immediate vicinity will provide key constraints on cosmological models of galaxy evolution (Springel et al. 2005;Croton et al. 2006).A crucial component in these numerical recipes is therefore the accurate determination of black hole masses as a function of cosmic history.
Robust black hole mass measurements have been limited mostly to nearby galaxies owing to the exquisite angular resolution needed to spatially resolve a black hole's sphere of influence.The technique of temporally resolving the structure of the broad-line region (BLR) around an actively accreting black hole, called reverberation mapping (RM; Blandford & McKee 1982;Peterson 1993), provides a viable alternative -or in the case of the distant universe, the only currently available -tool for directly measuring black hole masses.
Black hole scaling relations, such as that between the BLR size and the active galactic nuclei (AGN) restframe luminosity at 5100 Å (the radius-luminosity relationship; Wandel et al. 1999;Kaspi et al. 2000Kaspi et al. , 2005;;Bentz et al. 2006Bentz et al. , 2009aBentz et al. , 2013;;Shen 2013), enable "single-epoch" mass-determination methods for estimating black hole masses in broad-lined AGN out to high redshifts.The calibrations for such mass estimates, however, rest upon the presumptions that distant quasars are similar to the low-redshift reverberation-mapped AGN, while in fact they exhibit higher luminosities and larger Eddington ratios (Richards et al. 2011;Shen 2013).Whether virial assumptions for BLR gas dynamics based on local Seyferts can be extended to distant quasars must be tested against a broad range of black hole masses and luminosities.RM, as a tool that probes BLR gas structure, provides the means for such a test.
The RM method is successful for black hole mass determination because the emission lines in the BLR gas are found to respond to the stochastic continuum variations in an AGN via the transfer function (Peterson 1993) where L(v z , t) is the emission-line luminosity at line-ofsight velocity v z at observed time t, C(t) is the con-tinuum light curve, and Ψ(v z , τ ) is the transfer function that maps continuum variability to the emissionline response at v z after some time delay τ .This simple model assumes a linear response and no background light; fitting AGN light curves requires the linearized echo model that subtracts off the reference levels, e.g.L(v z , t) − L 0 (v z ) (see discussion in Horne et al. 2021).
Monitoring both the continuum and the emission-line light curves provides data that can allow the determination of Ψ(v z , τ ), also known as the velocity-delay map, whose shape depends on the structure and kinematics of the BLR (Horne et al. 2004).In practice, precise velocity-resolved information in the form of the transfer function demands intense monitoring with frequent sampling, high signal-to-noise ratio (S/N) data, and long duration.As a result, the number of AGN with such data available to-date has been limited.
With the goal of refining black hole scaling relations through expanding the local reverberation-mapped AGN database by which they are anchored, we embarked on a campaign aimed to monitor AGN at higher luminosity than samples targeted in the 2008 and 2011 campaigns of the Lick AGN Monitoring Project (LAMP, e.g., Bentz et al. 2009b; Barth et al. 2015).Higherluminosity AGN tend to have longer lags.Coupled with the fact that they have weaker variability amplitudes owing to the ample fuel supply in their accretion disks (vanden Berk et al. 2004;Wilhite et al. 2008;MacLeod et al. 2010), lag recovery is more challenging for these objects.The advent of large, multiobject programs such as SDSS-RM (Shen et al. 2015) and OzDES (Yuan et al. 2015) provides a new landscape where large numbers of reverberation lags for Hβ, Mg II, and C IV in quasars are determined over a broad redshift range (Shen et al. 2016;Grier et al. 2017Grier et al. , 2019)).Grier et al. (2017) measured Hβ time lags for 44 AGN with luminosities log[λL λ (5100 Å)/L ⊙ ] ≈ 43-45.5 at redshifts z = 0.12-1.Early results from the SDSS-RM campaign determined C IV lags in 52 quasars (Grier et al. 2019), whereas those from OzDES found two C IV-based black hole masses to be among the highest redshift (z = 1.9-2.6)and highest mass black holes [M BH = (3.3-4.4)×10 9 M ⊙ ] measured thus far with RM studies (Hoormann et al. 2019).Complementing these large multifiber spectroscopic studies, our campaign targeted AGN at low redshifts and aimed to yield high-fidelity data for velocity-resolved lag measurements and dynamical modeling.
During the past several years, different groups have determined resolved velocity-delay reverberation signatures for dozens of unique objects among local Seyferts, changing-look AGN, and those with Hβ asymmetry (e.g., Bentz et al. 2010a;Grier et al. 2013;Pancoast et al. 2014b;De Rosa et al. 2018;Du et al. 2018a;Williams et al. 2018Williams et al. , 2020;;Lu et al. 2019;Zhang et al. 2019;Sergeev 2020;Lu et al. 2021;Bentz et al. 2021).Our primary intent is to use a large, relatively broad sample to investigate statistical luminosity-dependent trends in BLR structure and gas dynamics using velocity-resolved reverberation, which could directly impact the accuracy of virial mass estimates.Here we present our first results that focus on new Hβ velocity-resolved measurements along with integrated Hβ lags and virial black hole masses.Our dataset will allow for forward-modeling work using, for instance, the CARAMEL code (Pancoast et al. 2011, 2014a, Villafaña et al., in prep.) to directly determine black hole masses.
This paper is organized as follows.The sample selection is described in Section 2. In Section 3, we present our observational program at Lick Observatory and details of the data reduction and processing work.Section 4 reports on our photometric campaign to provide the continuum light curves for our AGN sources.In Section 5, we illustrate our emission-line light curves and subsequent integrated and velocity-resolved Hβ lag detections, including an assessment of the lag significance and the observed variety of BLR kinematics.Section 6 depicts our line-width measurements and derived virial black hole masses, within the larger context of how our results compare with the existing AGN radiusluminosity relation.
We applied this 20-30-day criterion in Hβ lag or equivalently λL λ (5100 Å) to various catalogs of Type 1 AGN (e.g., Boroson & Green 1992;Marziani et al. 2003;Peterson et al. 2004;Vestergaard & Peterson 2006;Bachev et al. 2008;Winter et al. 2010;Shen et al. 2011;Zu et al. 2011;Joshi et al. 2012;Bennert et al. 2015;Sun & Shen 2015).Our selection was further narrowed with a redshift limit at z < 0.08 to ensure that the Hβ and [O III] lines fall blueward of the cutoff wavelength of the dichroic used in the Kast spectrograph (5500 Å).Sources were selected with declination δ ≥ −5 • and magnitude < 17 in the optical V or r band.We also assessed the short-term variability in previous light curves from the Catalina Real-time Transient Survey (Drake et al. 2009) for all our candidate sources and selected those with at least 0.1 magnitude of variations.Targets having high-quality velocity-resolved lag measurement from previous RM campaigns were excluded (e.g., Denney et al. 2010;Grier et al. 2012;Du et al. 2014;Wang et al. 2014), as were those no longer featuring broad hydrogen recombination lines in archival Lick optical spectra.
In order to ensure that lags can be measured accurately, we placed a further constraint requiring that the monitoring duration of an AGN, defined as the longest continuous period when the AGN can be observed at airmass < 2 during our campaign, be at least 3 times larger than the expected Hβ lag based on simulations.Given that our initial campaign duration was designed to be 9 months long (which extended to 1 year later on), this constraint resulted in a sample of 29 objects with duration-to-lag ratio > 3. Two sources with slightly shorter Hβ lag estimates (Mrk 315: 13 days; Mrk 704: 16 days) were included to better fill the right-ascension range of the sample.
After an initial probationary period of ∼ 1 month, we discarded 8 sources showing low continuum variability amplitude, retaining a sample of 21 objects for ongoing monitoring.This final sample of Seyfert 1s and their properties are listed in Table 1.

Overview of the Observational Program
We conducted a spectroscopic monitoring program of 1 yr duration at Lick Observatory on Mount Hamilton, California.Prior to the start of the campaign, we ran simulations to determine the minimum threshold for successful lag recovery among monitoring periods of varying lengths and cadences.Our estimated Hβ lag range required a sampling cadence of two nights per week over the course of a year to sufficiently resolve timedependent emission-line variations and to extend the temporal baseline for detecting robust reverberation sig-nal across the entire extent of the BLR, accounting for anticipated weather losses and seasonal gaps.
This project was allocated 100 nights at the Lick 3 m Shane telescope across three observing semesters between 2016 April 28 and 2017 May 6 (UT).The observing runs were distributed mostly with a frequency of ∼ 2 nights per week as requested, with the exception of bright time and occasional scheduling constraints owing to other time-sensitive observing programs.Partial nights were occasionally exchanged with or obtained from other programs to facilitate the sampling cadence of certain targets during this monitoring period.
Spectroscopic data were taken using the Kast Double Spectrograph (Miller & Stone 1994).
The red CCD detector was upgraded halfway into our campaign (September 17-20, 2016).The upgrade from the Reticon 400 × 1200 pixel CCD to the Hamamatsu 4096 × 1024 pixel CCD significantly improved the quantum efficiency, up to a factor of 2 at the red end, and removed severe fringing effects.Images taken with the new red CCD suffered from a high incidence rate of cosmic-ray hits and alpha-particle hits due to radioactive material in the new dewar window.These defects were removed at the reduction stage where frames were combined before the spectrum was extracted.
We employed the following setup for our AGN observations: the D55 dichroic with the 600/4310 grism (nominal coverage of 3300-5520 Å at 1.02 Å/pixel) on the blue side and the 600/7500 grating (nominal coverage of 4000-11,000 Å at 2.35 Å/pixel pre-upgrade, and 3800-10,000 Å at 1.31 Å/pixel post-uprade) on the red side, a compromise between broad spectral coverage and moderate spectral resolution.The plate scale was 0. ′′ 43/pixel for both the blue and the new red detectors, and 0. ′′ 78/pixel for the old red detector.A slit width of 4 ′′ was employed to minimize slit losses due to seeing variations between different nights (Filippenko 1982).A fixed position angle, optimized for each individual target, was chosen so that the spectroscopic aperture sampled the same portion of the host galaxy across the duration of the monitoring campaign.
Exposure times were typically up to 30 min per object each night, split into 2 or 3 separate exposures (preand post-upgrade of the red CCD, respectively) to facilitate cosmic-ray cleaning.Standard calibrations including bias exposures, arc frames (using the Ar, He, Hg, Cd, and Ne lamps), and dome flats were taken during the afternoon, and well-calibrated flux standard stars (G191B2B, Feige 34, BD+284211, HZ44, BD+262606, and/or BD+174708) were observed during the nights.Before the red CCD upgrade, red-side observations suffered from severe fringing at wavelengths longward of  7000 Å.In order to remove the fringe pattern, dome flats at the position of every object were taken before or after each standard star and AGN observation on the red side.This time-consuming procedure was no longer needed after the upgrade that eliminated the fringing effect.The observing parameters for our sample are listed in Table 2.
The wet 2016-2017 winter season at Lick hampered our observing effort, particularly during December through February.Overall, our spectroscopic campaign achieved an observational success rate of 70%, where 30 nights were lost entirely to weather or dome issues.For another 12 nights, only 6 or fewer AGN were observed when typically 12-15 objects would have been observed on a night with good conditions.We were confident that at least fourteen nights of observations were carried out under photometric conditions.Of the 70 nights when data were obtained, half were done with the old red CCD and half with the upgraded detector.
The number of total epochs observed for each object ranges from 22 to 50, with a median of 38 (see Table 2).

Spectroscopic Data Reduction
The Kast spectroscopic data were reduced using a combination of standard routines in IRAF and IDL, following the procedure outlined by Barth et al. (2015).
Here we provide a brief description of the reduction process of the blue-side data and pre-/post-upgrade red-side data.
In general, the data were processed via bias subtraction, flat fielding, cosmic-ray cleaning, one-dimensional (1D) extraction using an unweighted boxcar extraction region, wavelength calibration, and flux calibration.Error spectra were extracted and propagated through all subsequent calibration steps.Multiple spectra taken on the same night were combined.The width adopted for the spectral extraction was 10 ′′ , though it was widened to ∼ 15-20 ′′ for some sources (Mrk 315, Mrk 9, MCG Note-texp represents the typical total on-source exposure time for each AGN every night it was observed, usually split among 2-4 frames.N obs represents the total number of spectroscopic observations for each source.S/N represents the median S/N per pixel in the continuum at (5100-5200) (1 + z) Å. N phot represents the number of photometric nights for each source.
Most of the AGN were flux calibrated using the standard-star exposure that was closest in airmass.However, the flux-calibrated spectra for a few epochs exhibited some abnormal slopes and features close to the dichroic cutoff.We were unable to conclusively determine the cause of these anomalies, but they were most likely related to the dichroic.Since the [O III] λ5007 line, situated close to the dichroic cutoff, was assumed to be constant throughout the observing campaign and was used to normalize the nightly spectra, this peculiarity identified in several spectra added extra scatter to the Hβ light curves.In these cases, we attempted to correct the problem by calibrating the AGN with standard stars that were taken close in time for six nights.As a result, this flux-calibration issue was partially rectified but remained responsible for some residual scatter in the emission-line light curves For red-side data taken prior to 2016 Sep. 20, each observation was flattened using a dome-flat exposure taken at the same telescope position.For red-side data taken between 2016 Sep. 20 and 2017 Jan. 31, the Kast red dewar was slightly tilted and thus the two-dimensional spectra had to be adjusted with a rotation of 0.9 • before spectral extraction.

Photometric and Spectral Scaling
To ensure that all the spectra for each AGN are on consistent flux and wavelength scales across the monitoring period, a modified version of the procedure described by van Groningen & Wanders (1992, hereafter VGW92) was applied to the blue-side data for internal calibration (see Barth et al. 2015;Fausnaugh 2017, for more details).Having the spectra exhibit consistent spectral resolution throughout the temporal series is important particularly for extracting velocity-resolved measurements.First, we computed the shifts in wavelength among the time-series spectra from cross correlation.We applied these shifts to align the spectra and calculated an initial mean spectrum that was then designated as the reference.The [O III] λ5007 fluxes, presumed to be intrinsically constant throughout the course of the observing campaign, were measured from spectra taken during photometric nights and were averaged to estimate the true absolute flux of [O III].The number of photometric nights designated per object ranged from 3 to 12, with a median of 8 nights.The resulting reference spectrum was then scaled to match this measurement of the [O III] flux, with a median uncertainty of 8%.
Each night's spectrum was then aligned and scaled to the reference by matching the [O III] λ5007 emission-line profile following our modified VGW92 scaling method.In order to minimize variations outside of intrinsic AGN variability, we aimed to achieve a uniform spectral resolution across the full time series of spectra.Where VGW92 ignores resolution corrections for epochs observed at lower resolution than the reference spectrum, we adopted instead a modified approach similar to that described by Fausnaugh (2017): we first applied a Gaussian broadening kernel to the reference spectrum to ensure that all the rescaled spectra match a single resolution, typically the worst among the spectra.Excluding the occasional extreme outliers, this broadening kernel had a typical σ of ∼ 1.5-3 Å with a median σ of 2 Å for each AGN.Each of the spectra was then aligned by wavelength, scaled in flux, and broadened in spectral resolution to match the [O III] line profile of the broadened reference spectrum.
We performed a comparison of the modified VGW92 spectral scaling method to that using the Python code mapspec (Fausnaugh 2017) with the same wavelength windows and broadened reference spectrum.Differences between the two approaches include the following: (i) VGW92 calls for using a Gaussian kernel for smoothing, while mapspec offers an additional option of using Gauss-Hermite polynomials; and (ii) our modified VGW92 approach fits for free parameters via χ 2 minimization, while mapspec uses a Bayesian framework to optimize rescaling parameters and estimate model uncertainties.The runtime for the modified VGW92 method was typically ∼ 100 times shorter than for mapspec.We ran several tests comparing the level of scatter in the light curves of the [O III] line that resulted from both methods and concluded that modified VGW92 performed better for most of the AGN in our sample.We subsequently adopted the spectra scaled with the modified VGW92 approach for the ensuing analysis.
The scaling of the red-side data took a different approach.To normalize the red-side flux scale, red-side spectra were stitched to the corresponding blue-side spectra via aligning their overlapping regions (∼ 5300-5500 Å).All spectra were first aligned to the reference spectrum in wavelength.Because the ends of the spectra tended to be noisy, we used a weighted average to determine the overall multiplicative scaling factor.This technique worked well in most cases where the spectra were smooth or well behaved at the ends.The average scaled spectra for the sample are presented in Figure 1.

Noise-Corrected RMS Spectra
Relative variability across the spectrum can be visualized using the root-mean-square (rms) spectrum.The Peterson et al. (2004) procedure of taking the standard deviation of flux values at each wavelength element over all the epochs may inadvertently bias the rms spectrum given the inclusion of various noise factors in addition to genuine AGN variability (Barth et al. 2015).To remove the contribution due to photon-counting noise from the rms spectrum, Pei et al. (2017) suggested an approach to produce the "excess rms" (e-rms hereafter) spectrum defined per wavelength element λ and epoch i as where N is the number of epochs in the time series, F λ is the flux averaged over the time series, and F λ,i and δ λ,i refer to the flux and rms uncertainty in the flux at each λ and i, respectively.This was applied to the set of scaled spectra for each AGN.
In the case of an AGN with strongly varying broadline profiles, broad features are expected to dominate the e-rms spectra in the relevant spectral regions.The narrow [O III] lines, on the other hand, generally vanish in the rms given that they are assumed to be constant within the time-series data.Presented in Figure 2, the e-rms spectra generally feature a strong blue continuum (except for Mrk 315), broad-line emission in the Balmer series and (usually) in He Here, as in Table 1, the AGN are shown in order of increasing right ascensions.
residual narrow [O III] features confirms that our internal flux calibration worked well near this wavelength, but other residual narrow-line features are likely artifacts of imperfect flux calibration errors that increase at bluer wavelengths.This is particularly the case for Mrk 315, which has a red e-rms continuum slope and prominent [O II] λ3727.

Spectral Decomposition
The traditional approach of measuring broad emission-line fluxes using a simple linear continuum subtraction is subject to several inadequacies.A linear continuum model is unable to separate blended emissionline features, such as He II or Fe II lines that often overlap the broad Hβ profile.Furthermore, a linear fit is a poor approximation to the actual continuum, which includes both AGN and host-galaxy components.Residual errors from continuum subtraction can be particu-larly problematic for velocity-resolved RM in the faint high-velocity wings of broad emission lines, where oversubtraction or undersubtraction of the continuum could lead to biased inferences on BLR structure and kinematics.In recent years, spectral decomposition approaches have been applied to the data from some RM campaigns (e.g., Barth et al. 2013Barth et al. , 2015;;Hu et al. 2015).By fitting a multicomponent model to each night's spectrum, the contributions of individual line and continuum components can be better isolated.We find that in most of our objects, the Hβ line profiles isolated using the two methods do not differ dramatically except in the red wings, but spectral decomposition is useful for separating blended line components and removes the starlight component for objects with high starlight fraction much more effectively.
Here we have adopted and applied the spectral fitting method described by Barth et al. (2015)  Hβ spectral region.We refer the reader to Barth et al. (2015) for a detailed description, and provide a brief overview of the method here.The model components include a power-law AGN continuum, starlight from a single-burst old stellar population at solar metallicity (Bruzual & Charlot 2003) convolved with a Gaussian velocity broadening, and emission lines including [O III] (narrow), Hβ (broad and narrow), He II (broad and narrow), He I (broad), and an Fe II emission template convolved with a Gaussian velocity broadening.Emission-line profiles were modeled using a fourth-order Gauss-Hermite function, except for the He I and He II lines for which a Gaussian was used.Each line's velocity centroid was allowed to vary independently, except for the [O III] λλ4959, 5007 lines which were required to have the same velocity profile and a 1:3 flux ratio.For the Fe II blends, we tested different template spectra as described by Boroson & Green (1992), Véron-Cetty et al. (2004), and Kovačević et al. (2010).In the near future, we will employ a promising, new Fe II template spectrum based on Mrk 493 (Park et al., in prep.).Similar to Barth et al. (2015), the multicomponent Kovačević et al. (2010) template provided the best fit to the Fe II lines in all of our AGN and was used for the final fits, except for the case of I Zw 1.For I Zw 1 we found that the Boroson & Green (1992) template, which is based on I Zw 1 itself, provided the best fit with minimized residuals.As part of the fitting process, a Cardelli et al. (1989) reddening model is applied to the model spectrum, allowing E(B − V ) to be a free parameter.Fitting was carried out over a rest-wavelength range of approximately 4200- 5200 Å, with the exact range tailored to the data for each AGN depending on its redshift.The Hγ + [O III] λ4363 blend was masked out from the fit, because decomposing this blend would add several additional free parameters to the model.The full model, including the multicomponent Kovačević et al. (2010) iron template, includes 33 free parameters.Model fits were optimized using a Levenberg-Marquart algorithm as implemented by Markwardt (2009).For each AGN, the model was first fitted to the mean spectrum, and the fit parameters from the mean spectrum fit were used as the starting parameter estimates for the fit to each individual night's spectrum.The overall mean spectrum as fitted with the different model components for each galaxy is illustrated in Figure 3.

THE PHOTOMETRIC CAMPAIGN
AGN continuum light curves were measured from imaging data.We chose the V band since it provides a fairly clean continuum measurement with relatively little contamination from broad emission lines for lowredshift sources.To maximize the chances of detecting the delayed response of the emission-line variations to those of the continuum, we began the photometric monitoring campaign two months before the start of the spectroscopic program.From February 2016 to May 2017, our team used a network of eight telescopes across the world and obtained high-fidelity V -band images with up to nightly cadence within each object's monitoring season.

Telescopes and Cameras
Our photometric campaign included observations with the following telescopes: the 0.76 m Katzman Automatic Imaging Telescope (KAIT) (Filippenko et al. 2001) and the Anna Nickel telescope at Lick Observatory on Mount Hamilton, California; the Las Cumbres Observatories Global Telescope (LCOGT) network (Brown et al. 2013;Boroson et al. 2014); the 2 m Liverpool Telescope at the Observatorio del Roque de Los Muchachos on the Canary island of La Palma, Spain (Steele et al. 2004); the 1 m Illinois Telescope at Mount Laguna Observatory (MLO; Smith & Nelson 1969) in the Laguna Mountains, California; the San Pedro Mártir Observatory (SPM) 1.5 m Johnson telescope (Butler et al. 2012;Watson et al. 2012) at the Observatory Astronómico Nacional located in Baja California, México; the Fred Lawrence Whipple Observatory 1.2 m telescope on Mount Hopkins, Arizona; and the 0.9 m West Mountain Observatory (WMO) telescope at Utah Lake in Utah.KAIT, LCOGT, SPM-1.5 m, and Liverpool are fully robotic telescopes.Table 3 summarizes the telescope and detector properties.

Photometry Measurements and Continuum Light
Curves All images were processed with reduction procedures applied as part of the standard pipeline for each facility, including overscan subtraction and flat-fielding.For telescopes that did not include world coordinate system (WCS) solutions as part of their standard processing, we used astrometry.net(Lang et al. 2010) to add WCS information to the FITS headers.
Measurement of the AGN light curves was carried out using the automated aperture photometry pipeline described by Pei et al. (2014), which can be used to measure AGN light curves using data from any number of telescopes with diverse camera properties.This procedure, written in IDL and based on the photometry routines in the IDL Astronomy User's Library (Landsman 1993), is designed to automate the process of identifying the AGN and a set of comparison stars in each image by their coordinates, measuring their instrumental magnitudes, and using the comparison-star measurements to obtain a consistent magnitude scale across the full time series for data from multiple telescopes.
The aperture-photometry routine returns photometric uncertainties (in magnitudes) based on the photoncounting statistics, background uncertainty, and CCD readout noise.However, other error sources contribute to the actual error budget, such as imperfect flat-fielding and point-spread function (PSF) variations across the field of view.To account for these additional error sources, we measured the excess variance σ 2 x in the light curves of the comparison stars for each AGN, for each telescope's data.Assuming that the excess variance in the comparison-star light curves is a good estimate of the additional error over and above the statistical uncertainties, we inflated the fractional flux uncertainties on the AGN photometry by addition in quadrature using σ 2 tot = σ 2 phot + σ 2 x , where σ phot is the original photometric uncertainty on a given data point, and σ 2 x is the average (over all comparison stars in the field) of the excess variance in the comparison-star light curves.This error adjustment was carried out separately for each telescope's data.
AGN light curves measured from different telescopes tend to have small offsets in magnitude relative to one another as a result of differences in filter transmission curves, even after normalizing them to the same set of comparison stars.To correct for these offsets, we applied an additive shift (in magnitudes) to bring each telescope's data into best average agreement with the Rest Wavelength (Å)  data from the telescope having the most data points for each AGN.Finally, the V -band magnitude scale for each AGN field was calibrated using observations of standard stars taken on photometric nights at WMO and crosschecked against the AAVSO Photometric All-Sky Survey (Henden et al. 2016).
Magnitudes were converted to f λ for measurement of broad emission-line lags.The typical uncertainties of the photometric data points δ f are at the 0.4% level, not exceeding 0.7% for any AGN.The final V -band light curves are presented in Table 4 and shown alongside the emission-line light curves in §5.

EMISSION-LINE LIGHT CURVES AND LAG DETERMINATION
We used results from spectral decomposition to derive Hβ emission-line light curves.The flexibility of the multicomponent modeling afforded us a number of ways to extract the Hβ flux.In the spectral region near the Hβ line, sources with strong Fe II emission were subjected to substantial degeneracy of spectral decomposition that could result in, from night to night, instability in model fitting in terms of how much flux gets assigned to the different model components.This instability results in additional noise in the Hβ light curve.
To minimize the noise due to degeneracy from spectral fits, we used a version of the decomposed model spectra where only the AGN and stellar continuum components have been removed.This version consisted of the residuals of all present emission lines after continuum sub-traction.We then selected a corresponding extraction window free of other emission lines for the calculation of Hβ flux (see Table 5).A different way to extract Hβ flux is to use a version of the decomposed model spectra that only contained Hβ broad and narrow components, with other emission-line components as well as the AGN and stellar continua removed.The resulting Hβ light curves computed from the two methods were very similar for objects that have relatively clean Hβ spectral regions, but those for the emission-line residual versions were less noisy for objects with more complex line profiles within the Hβ spectral region.The final Hβ light-curve measurements for the full sample are presented in Table 6.The light curves for the V band and Hβ are shown in Figures 4-10.The rest of the emission-line features will be presented in a forthcoming paper.To account for the residual uncertainties from spectral scaling, we computed the excess variance in the [O III] light curve for each AGN as the additional uncertainty term.We then combined this excess variance term with the statistical error in quadrature such that Both the statistical and the modified errors are listed in Table 6.Following Bentz et al. (2009a), we compute the variability statistics F var for each light curve and found that it ranges from 0.028 to 0.268, with a median of 0.069 (see

Integrated Hβ Lags
Following White & Peterson (1994), Peterson et al. (2004) and others, cross-correlation was employed to determine the delay in the integrated Hβ line signal relative to that from the V -band continuum from the photometric campaign.We chose to adopt the interpolated cross-correlation function (ICCF; e.g., Peterson et al. 1998) method.Given a lag range that represents reasonable bounds of the expected lag, the cross-correlation function (CCF) between the two light curves is com-puted at every time lag τ at 0.5 day intervals via ) where L ( C) and σ L (σ C ) are the mean and standard deviation of the emission-line (continuum) time series, respectively.
Two lag measures are computed: τ peak , the lag at the peak of the CCF, and τ cen , the centroid of the CCF for all points in the CCF above a predetermined threshold value (80% of the peak CCF value).We adopt the latter for the subsequent analysis.For error analysis, we employed the Monte Carlo flux randomization method (Peterson et al. 1998) and repeated the CCF calculation for 10 3 realizations, building up distributions of correlation measurements.The median and ±1σ widths of the cross-correlation centroid distributions (CCCD) were then adopted as the final lag measurements, and their associated uncertainties as shown in Figures 4-10.

Assessment of Cross-Correlation Reliability
In any RM campaign, it is expected that some AGN will yield highly robust measurements of reverberation lag, thanks to strong variability and high-quality data, while other objects may exhibit little or no evidence for correlated variability between the continuum and emission-line light curves.This can occur as a result of observational factors including low S/N or poor temporal sampling, or factors intrinsic to the AGN including low variability amplitude or low responsivity of the emission line to variations in the ionizing continuum.There is not necessarily a clear demarcation between reliable and unreliable lag measurements, and a variety of methods have been used to assess the significance or quality of lag detections (e.g., Grier et al. 2017).In some cases, a simple threshold value of the correlation strength r max is used, but r max alone is not necessarily a good indicator of correlation significance for red-noise light curves.An alternative method is to employ null-hypothesis testing to determine the probability that two uncorrelated red-noise light curves having the same S/N and cadence as the data would yield an r max at least as strong as the observed value.Such methods have been employed for analysis of multiwavelength continuum correlations in AGN (e.g., Uttley et al. 2003;Arévalo et al. 2008;Chatterjee et al. 2008), but until recently have rarely been employed to assess broad-line reverberation lags (Penton et al. 2021;Li et al. 2021).Our method will be presented in detail by Guo et al. (in prep.), and we provide a brief description here.U et al.We first generated light curves according to the damped random walk (DRW) model with the Python software CARMA1 (Kelly et al. 2009(Kelly et al. , 2014)).
The DRW model provides an adequate description of ultraviolet/optical variabilities in AGN, though with plausible deviations on short (McHardy et al. 2006;Mushotzky et al. 2011;Kasliwal et al. 2015;Smith et al. 2018) or long (MacLeod et al. 2010;Guo et al. 2017) timescales.Each simulated light curve represents a segment randomly selected from a 100-times longer light curve predicted from the same DRW model fitted to the observed light curve.We resampled each mock light curve to have the same cadence as the real observations and added Gaussian random noise based on the S/N of the data, creating 10 3 realizations of each light curve.Cross-correlation measurements were then carried out using the simulated data to determine the distribution of r max values that would be obtained for uncorrelated light curves.A two-way simulation was performed -that is, we calculated the CCF between the real continuum and each simulated emission-line light curve for the first 500 simulations, and then the opposite way for the rest of the realizations.
We measured the CCF for all the simulations, using the same lag search range employed for each AGN, and counted, out of 10 3 , the number of positive lags (τ > 0) with peak values r max higher than our observed r max .The resulting fraction represents the probability that a correlation signal found between two uncorrelated light curves would exceed the correlation signal of the data.This derived p-value, denoted p(r max ), thus provides an indicator of the robustness of our lag detections given the observed r max and the observed properties of the light curves including S/N, cadence, and duration.
We emphasize that the p-values derived from this method are not false-positive probabilities, and they do not give a measure of the absolute "significance" of a lag detection.Strictly speaking, our p-values simply give the probability that uncorrelated light curves having the same statistical properties as the data would give a cross-correlation signal as strong as that seen in the data.Consequently, a smaller p-value denotes a more robust and reliable detection of the correlation signal between two light curves.In general, for broad-line RM we have a strong prior that a reverberation signal is very likely to be present in the data, thus our null hypothesis of intrinsically uncorrelated light curves is an extreme and fairly unlikely scenario.We further emphasize that there is no strict cutoff between significant and insignif-icant lag detections by this method; the p-value merely gives an indication of the relative degree of reliability between different measurements.
We assess the quality of our resulting lags based on these quantitative assessment indicators in Figure 11.As expected, there is an anticorrelation between r max and p(r max ), but it is not a tight one-to-one relationship owing to the differences in S/N, sampling cadence, and variability amplitude among our sample.The results are consistent with the general expectation that the highest-quality light curves that exhibit clear correlated variability have high r max and small p(r max ) values.Considering these lag-assessment results, we conclude that 16 of the 21 AGN in our sample falling in the region p(r max ) < 0.2 and r max > 0.6 have correlations between continuum and Hβ light curves that are sufficiently robust for the ensuing analyses (noting that they display a broad range of p-values and thus range from highly robust to relatively weak detections of correlated variability), while we discard the remaining five objects from further analysis based on their low r max and high p-values.Based partly on visual inspection of the light-curve quality, our chosen r max threshold is relatively conservative compare to those selected by other large surveys (e.g., minimum r max = 0.45 by SDSS-RM; Grier et al. 2017).These r max and p(r max ) values, along with AGN luminosities and lag measurements τ cen in both observed and rest frames, are reported in Table 7.

Velocity-Resolved Hβ Lags
To obtain velocity-resolved reverberation results as a means to probe the kinematics of the BLR gas, we used CCF to compute lag measurements for individual velocity segments of the emission line.Our procedure follows that described by Bentz et al. (2009b), Denney et al. (2009), andGrier et al. (2013), and is detailed below.
To verify the effect of binning on the apparent BLR kinematics, we divided the broad Hβ component into eight bins via two different schemes: (i) bins of equal rms flux, and (ii) bins of equal velocity width.In the first approach, we initially determined the total rms flux from integrating the Hβ line in the e-rms spectrum.Then we established velocity bins such that the rms flux within each bin would equal one-eighth of the total flux.There are some exceptions in the cases where the red wing of the rms profile suffered from significant noise and the last bin was thus truncated at the red line wing.In the second approach, bins are divided evenly across the width of the line in velocity space.In either binning scheme, light curves were computed from each bin of the continuum-subtracted Hβ profile and lags were measured against the V -band continuum light curve.Resulting lags from the two binning schemes generally agree in the line's core where the S/N is high; minor disagreement is found in the noisy line wings of some sources.For clarity, we adopt the first scheme (equal-rms-flux binning) and show the resulting velocity-resolved lag spectra for each AGN with robust lag in Figure 12 and in Table 8.

Implications for BLR Kinematics
The velocity-resolved lags allow us to make simple qualitative inferences about the BLR geometry and kinematics from the transfer-function models (e.g., Welsh & Horne 1991;Horne et al. 2004;Bentz et al. 2009b).
Qualitatively, symmetric velocity-resolved structure around zero velocity is consistent with either Keplerian, disk-like rotation or random motion without net radial inflow or outflow over an extended BLR.In such a case, the high-velocity wings exhibit the shortest lags because the highest rotation speeds correspond to material orbiting closest to the central black   hole.Radially infalling gas produces an asymmetric pattern displaying longer lags at the high-velocity blue wing, while the opposite is true in the case of outflowdominated kinematics.We note that the presence of absorption, noise, or other geometric complexities may muddle these interpretations, which could be elucidated with direct modeling using forward-modeling codes such as CARAMEL in follow-up work (Villafaña et al., in prep.).We find a diversity of velocity-resolved lag spectra among our sample representing the qualitative signature of symmetric (3C 382, MCG +04−22−042, Mrk 110, Mrk 1392), infalling (Mrk 1048, Mrk 704, Zw 535−012, Mrk 841, RXJ 2044.0+2833,NPM1G+27.0587,SBS 1518+593), and outflowing (RBS 1303 andRBS 1917) scenarios.In the cases of Ark 120, Mrk 9, and PG 2209+184, the velocity-resolved structure appears flat across the Hβ profile within the uncertainties, which may be difficult to describe with simple models.Whether these kinematics exhibit luminosity-dependent trends and how a subset of these results relate to other RM samples will be further examined in §6 and the Appendix, respectively.

AGN Radius-Luminosity Relationship
The applicability of the radius-luminosity relationship to single-epoch mass determination (Shen et al. 2008;Shen & Liu 2012) relies on the small scatter initially quantified (∼ 0.2 dex; Kaspi et al. 2005;Bentz et al. 2009aBentz et al. , 2013;;Kilerci Eser et al. 2015;Martínez-Aldama et al. 2020), but recent studies of in-creasingly diverse AGN samples led to an increase in this scatter (e.g., Grier et al. 2017;Du et al. 2018b).The main reason for this increased scatter seems to depend on the Eddington ratio (Bian et al. 2012;Du et al. 2015Du et al. , 2016a;;Martínez-Aldama et al. 2019;Dalla Bontà et al. 2020), an important third parameter suggesting that there is more diversity in the AGN population than can be explained by the simple two-parameter radiusluminosity relationship.
Here we compare our sample on the AGN radius-luminosity relationship with other RM samples from the literature (Bentz et al. 2013;Grier et al. 2017;Du et al. 2018b) in Figure 13.The best-fit line from SDSS-RM for this relation is log(R BLR /ltday) = K + α log(λL λ (5100 )/10 44 erg s −1 ), (5) with a median slope α = 0.45 and a median normalization K = 1.46 (Fonseca Alvarez et al. 2019).We computed the AGN luminosity λL λ (5100 Å) based on the AGN continuum flux at 5100 Å as extracted from the decomposition of the mean spectrum.Luminosity distances were derived using the cosmology calculator (H 0 = 67.8km s −1 Mpc −1 ) provided by Wright (2006).All of our AGN fall within 0.5 dex of the radiusluminosity relation as defined by the literature points within the range of λL λ (5100 Å) = 10 43.5−44.4erg s −1 .The best-fit line in the form of Equation 5 for all the data points combining our current work with those from the literature has a normalization K = 1.35 and α = 0.38.
We further explore whether BLR kinematics might depend on AGN luminosity.Past work has shown that while inflow or outflow are indicated in a subset of AGN, the majority of Hβ velocity-delay maps appear roughly symmetric with rotationally dominated kinematics (e.g., Bentz et al. 2009bBentz et al. , 2010b;;Denney et al. 2010;Barth et al. 2011;Horne et al. 2021).The AGN with velocity-resolved reverberation detected to-date are mostly within the low-luminosity range [λL λ (5100 Å) = 10 42−43.5 erg s −1 ].Our results add significantly to the existing number of sources having velocity-resolved measurements with λL λ (5100 Å) ≈ 10 44 erg s −1 .While we caution that velocity-resolved lag spectra provide only a qualitative impression of the BLR kinematics, we illustrate the BLR kinematics via different colored symbols shown alongside those from the literature in Figure 14.No trend is seen between BLR kinematics and AGN luminosity either within our LAMP2016 sample or when combined with past results from the literature.It is plausible that the BLR shows timedependent structural changes, from both the kinematic and the geometric points of view, and there may be a lag in its correlation with the luminosity of the AGN.The diversity in BLR kinematics found in our mod-  Bentz et al. 2013;Grier et al. 2017;Du et al. 2018b).The best-fit relation from SDSS with a median slope α of 0.45 and a median normalization K of 1.46 (Fonseca Alvarez et al. 2019) is plotted as the dashed line.The best-fit line in the form of Equation 5 for all the data points combining our current work with those from the literature (blue) has a normalization K = 1.35 and α = 0.38.
erately high-luminosity sample appears consistent with that seen among AGN with super-Eddington accretion rates (e.g., Du et al. 2016b).

Hβ Line Widths
Following the procedure outlined by Barth et al. (2015), we measured the line width of the broad Hβ component from both the spectrally-decomposed mean spectrum and the rms spectrum.The mean spectrum constructed from the time-series data has a continuum level of zero with just residual noise from spectral fitting, but the continuum of the rms spectrum incorporates photoncounting errors in addition to residual noise.The noise continuum is removed after a simple linear fit before making line-width measurements.
We computed the line width using two separate parameters according to Peterson et al. (2004): the full width at half-maximum intensity (FWHM) of the line, and the line dispersion σ line .The FWHM of a singlepeaked line is defined as the difference between the wavelengths blueward and redward of the peak P (λ) max that correspond to half of its height.A double-peaked line requires additional attention to defining P (λ blue,red ) max for the blue and red peaks, respectively, and thus their corresponding wavelengths at half the maximum height, but otherwise the procedure in determining the differ- ence between the resulting wavelengths is similar.As for σ line (square root of the second moment of the Hβ line profile), we adopted from Peterson et al. (2004) the following: where and In each of these calculations, we also varied the endpoints of the Hβ extraction window randomly within ±5 Å ten times in order to account for uncertainties brought about by the precise choice of the Hβ spectral limits.The resulting line-width measurement is the mean from these ten trials.
For error analysis associated with the measured width parameters, we undertook a Monte Carlo bootstrapping procedure and simulated 100 realizations of the mean and rms spectra, respectively, based on random subset sampling of the existing nightly spectra.For a source with n epochs of observations, we simulate new mean and rms spectra of the broad Hβ component from combining n randomly-drawn time-series spectra with repetition allowed.We repeated our measurement of the FWHM and σ line on these new spectra, where the means and the standard deviations of the resulting distributions were taken as the line widths and their corresponding uncertainties, respectively.
The final line-width measurements are corrected for instrumental broadening subtracted in quadrature from the observed line width ∆λ obs : Given that we used the same blue-side spectral setup as the LAMP2011 campaign, we adopt from Barth et al. (2015) the instrumental FWHM of 380 km s −1 and dispersion σ inst ≈ 162 km s −1 for our small corrections.
The resulting FWHM and σ line measured from both the mean and the e-rms spectra are presented in Table 9.

Virial Black Hole Masses
The mass of a black hole may be determined most robustly using dynamical modeling codes such as CARAMEL (e.g., Pancoast et al. 2011).Nonetheless, virial estimates are useful within uncertainties that depend on the BLR geometry.The black hole mass and its associated uncertainty are computed and propagated from the virial equation where f is the scaling factor that depends on the geometry and kinematics of the BLR, c is the speed of light, τ ± σ τ is the time delay with uncertainties, v ± σ v is the velocity of the BLR gas with uncertainties as measured by the line width, and G is the gravitational constant.The empirically fitted virial factor from Woo et al. (2015) is consistent with that from CARAMEL dynamical modeling (Pancoast et al. 2014b), so we have adopted for our black hole mass calculations f = 10 0.65 = 4.47 from Woo et al. (2015), τ = τ cen , and v = σ line (rms) from Dalla Bontà et al. ( 2020) to facilitate comparisons with other studies.We present both the virial products cτ cen σ 2 line /G and the derived f -dependent black hole masses in Table 10.Comparisons of our results to existing prior black hole mass measurements for several sources can be found in the Appendix.Overall, we find agreements between our measurements and those from the literature within uncertainties.

CONCLUSIONS
We present the first results from our 100-night LAMP2016 campaign where we monitored a sample of 21 luminous Seyfert 1 nuclei with AGN luminosity λL λ (5100 Å) ≈ 10 44 erg s −1 during April 2016 -May 2017 at Lick Observatory.Our analysis here provided the Hβ emission-line light curves and integrated Hβ lag detections measured against V -band continuum photometric light curves for the full sample.We further assessed the significance of the lag determinations, and computed, for a subset of sources with goodquality lags, their velocity-resolved reverberations, inferred BLR kinematics, broad Hβ line widths, and virial black hole mass estimates.
The unusually rainy weather at Mount Hamilton during the winter months of 2016-2017 hampered our ability to monitor the variability of our AGN as closely as would have been ideal, leaving large gaps in several of our emission-line light curves.Consequently, we were not able to measure reliable lags for six of our AGN owing to inadequate light curves.Our overall results nearly double the number of existing velocity-resolved lag measurements, particularly at the moderately high-luminosity regime among previous reverberation-mapped samples, and revealed a diversity of signatures potentially indicative of different BLR gas kinematics.Given the complexity of the multiple factors involved, the lack of any clear correlation between AGN luminosity and velocity-resolved lag structure suggests that luminosity itself is not the sole factor responsible for setting the kinematic state of the BLR.It may also depend on other properties of the accreting black holes.Follow-up direct dynamical modeling work will shed light on the detailed kinematics and provide robust dynamical black hole masses to help further calibrate widely-adopted black hole scaling relations.This work makes use of observations from the LCOGT network.The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council.Based on observations acquired at the Observatorio Astronómico Nacional in the Sierra San Pedro Mártir (OAN-SPM), Baja California, México, we thank the daytime and night support staff at the OAN-SPM for facilitating and helping obtain our observations.Some of the data used in this paper were acquired with the RATIR instrument, funded by the University of California and NASA Goddard Space Flight Center, and the 1.5-meter Harold L. Johnson telescope at the Observatorio Astronómico Nacional on the Sierra de San Pedro Mártir, operated and maintained by the Observatorio Astronómico Nacional and the Instituto de Astronomía of the Universidad Nacional Autónoma de México.Operations are partially funded by the Universidad Nacional Autónoma de México (DGAPA/PAPIIT IG100414, IT102715, AG100317, IN109418, IG100820, and IN105921).We acknowledge the contribution of Leonid Georgiev and Neil Gehrels to the development of RATIR.This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF grant AST-1412587 and contributed by observers worldwide.We acknowledge the use of The AGN Black Hole Mass Database as a compilation of some of the reverberation mapped black hole masses prior to 2015 (Bentz & Katz 2015).This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Figure 1 .
Figure1.Mean spectra for the LAMP2016 sample including both blue-and red-side data.Here, as in Table1, the AGN are shown in order of increasing right ascensions.

Figure 2 .
Figure 2. The blue-side e-rms spectra for the VGW92-scaled time-series data featuring strong blue continuum and broad-line emission in the Balmer series, overlaid with dashed lines indicating the wavelengths of the prominent emission features (labeled in the top panels).The absence of residual narrow [O III] features demonstrates the effectiveness of our internal flux calibration at the red end, but other residual narrow-line features (such as the [O II] λ3727 feature in Mrk 315) is most likely due to a combination of underlying continuum variability and imperfect spectrophotometric calibration errors at the blue end.

Figure 3 .
Figure 3. Spectral decomposition results.The observed mean spectrum (black) for each galaxy is plotted alongside the decomposed model components: starlight (purple), AGN power-law continuum (green), Hβ (magenta), He II (cyan), He I (orange), Fe II (grey), and [O III] (blue).The sum of the fits is represented by the red line.The He I component is negligible in most cases.The blend of Hγ and [O III] λ4363 is excluded from the fit.The Kovačević et al. (2010) Fe II template was applied to all sources except for I Zw 1, which was fitted with the Boroson & Green (1992) template.

Figure 4 .
Figure 4. Light curves for the LAMP2016 sample (Zw 535−012, I Zw 1, Mrk 1048): V -band continuum flux density f λ in 10 −15 erg cm −2 s −1 Å −1 , color-coded by telescopes (top left); the Hβ emission line flux f line in 10 −15 erg cm −2 s −1 (bottom left).The error bars plotted here for the Hβ light curves incorporate the uncertainty term from the normalized excess scatter of [O III].The cross-correlation function is shown in the right panel for each AGN, alongside the cross-correlation centroid distribution in yellow.The dotted vertical line indicates the median value of the CCF.

Figure 11
Figure 11.p(rmax) versus rmax as two objective lagassessment parameters for our sample.Sources with p(rmax) above 0.2 or rmax less than 0.6 are deemed to have unreliable lags as represented by open circles.

Figure 12 .
Figure 12.Velocity-resolved reverberation lags for our sample.In each of the top panels, the Hβ lags in days (y-axis) shown are for bins of equal rms flux, plotted against relative velocity in km s −1 (x-axis).The blue dotted line indicates the overall integrated Hβ lag with ±1σ uncertainty spanned by the blue bar.The middle and bottom panels show the Hβ profile of the mean and rms spectra, respectively.

Figure 13 .
Figure 13.Radius-luminosity relationship of our sample (blue, filled) superimposed on top of other work from the literature (open symbols;Bentz et al. 2013;Grier et al. 2017;Du et al. 2018b).The best-fit relation from SDSS with a median slope α of 0.45 and a median normalization K of1.46 (Fonseca Alvarez et al. 2019) is plotted as the dashed line.The best-fit line in the form of Equation 5 for all the data points combining our current work with those from the literature (blue) has a normalization K = 1.35 and α = 0.38.

Figure 14 .
Figure 14.Radius-luminosity relationship of our sample (filled and color-coded by BLR kinematics) along with other velocity-resolved results from Bentz et al. (2010a), Pancoast et al. (2014b), Grier et al. (2013), De Rosa et al. (2018), and Du et al. (2018b) (open symbols with the same color coding).No trend is seen for the different BLR kinematics within our LAMP2016 sample or when combined with the literature.
. The Kast red CCD detector upgrade, led by B. Holden, was made possible by the Heising-Simons Foundation, William and Marina Kast, and the University of California Observatories.KAIT and its ongoing operation were made possible by donations from Sun Microsystems, Inc., the Hewlett-Packard Company, AutoScope Corporation, Lick Observatory, the NSF, the University of California, the Sylvia & Jim Katzman Foundation, and the TABASGO Foundation.Data presented herein were obtained using the UCI Re-mote Observing Facility, made possible by a generous gift from John and Ruth Ann Evans.Research at UC Irvine has been supported by NSF grants AST-1412693 and AST-1907208.V.U acknowledges funding support from the University of California Riverside's Chancellor's Postdoctoral Fellowship and NASA Astrophysics Data Analysis Program Grant #80NSSC20K0450.Her work was conducted in part at the Aspen Center for Physics, which is supported by NSF grant PHY-1607611; she thanks the Center for its hospitality during the "Astrophysics of Massive Black Holes Merger" workshop in June and July 2018.T.T. acknowledges support by the Packard Foundation through a Packard research fellowship.V.N.B. and I.S. gratefully acknowledge assistance from NSF Research at Undergraduate Institutions (RUI) grants AST-1312296 and AST-1909297.Note that findings and conclusions do not necessarily represent views of the NSF.G.C. acknowledges NSF support under grant AST-1817233.J.H.W. acknowledges the funding from the Basic Science Research Program through the National Research Foundation of Korean Government (NRF-2021R1A2C3008486).

APPENDIXA .
COMPARISON WITH PRIOR M BH MEASUREMENTS RM studies over the years have shown the intrinsic scatter among BH mass measurements from the RM technique to be inherently small, demonstrating its robustness for measuring BH masses.In the well-studied case of NGC 5548 with many epochs of RM-based BH masses (see Dalla Bontà et al. 2020, for a recent compilation), the distribution of U et al. 42.1 erg s −1 and FWHM(Hβ) = 2555 km s −1 from a single-epoch spectroscopic observation in 2006.

Table 1 .
Sample Properties Note-Objects in this and subsequent tables and figures are listed in RA order.Redshifts are from NED. References for λL λ (5100 Å) used to estimate the Hβ lag in our sample selection: (1) Previous Lick spectra; (2) Marziani et al.

Table 3 .
Telescope and CCD Properties for Imaging Observations

Table 4 .
V -band Light Curve Data Note-Dates are listed as HJD−2450000.Units for continuum flux density f λ are 10 −15 erg s −1 cm −2 Å −1 .This table is published in its entirety in machine-readable format on the online version of the journal.Sample entries are shown here for guidance regarding the table's form and content.

Table 6 .
Hβ Light Curve Data Note-Dates are listed as HJD−2450000.Units for emission-line flux f are in 10 −15 erg s −1 cm −2.The original δ f represents the photon-counting errors on the flux measurement, while the modified δ f incorporates the spectral scaling uncertainties from [O III] scatter.This table is published in its entirety in machine-readable format on the online version of the journal.Sample entries are shown here for guidance regarding the table's form and content.

Table 7 .
AGN Luminosity and Hβ Cross-Correlation Lag Results

Table 9 .
Rest-Frame Line Widths of the Broad Hβ Component Note-All line-width measurements are in km s −1 .theLickobserving campaign: (Shane) Zachary Parsons, Estefania Padilla Gonzalez, Noah Rivera, Cristilyn Gardner, Jake Haslemann, Sean Lewis, and Ellen Glad; (Nickel) Nick Choksi, Sameen Yunus, Jeff Molloy, Andrew Rikhter, and Haynes Stephens.Photometric data collection at MLO was supported by NSF grant AST-1210311; we thank Robert Quimby, Emma Lee, Joseph Tinglof, Eric McLaughlin, Amy Igarashi, and Tariq Johnson for assistance with these observations.We are deeply grateful to the UCO/Lick staff for help with scheduling and supporting the observations.Research at Lick Observatory is partially supported by a generous gift from Google

Table 10 .
(Woo et al. 2015and Derived BH MassesNote- † Assuming f = 4.47(Woo et al. 2015).A.V.F.'s group at U.C. Berkeley is grateful for support from the TABASGO Foundation, the Christopher R. Redlich Fund, the Miller Institute for Basic Research in Science (in which he is a Miller Senior Fellow), and many individual donors.K.H. acknowledges support from STFC grant ST/R000824/1.We acknowledge the generous support of Marc J. Staley, whose fellowship partly funded B.E.S. whilst contributing to the work presented herein as a graduate student.I.S. acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306.Research by S.V. is supported by NSF grants AST-1813176 and AST-2008108.