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

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 1 yr to sufficiently resolve time-dependent emission-line variations and to extend the temporal baseline for detecting a robust reverberation signal 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 (2016 September 17–20). 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 two 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–11000 Å at 2.35 Åpixel⁻¹ pre-upgrade, and 3800–10000 Å at 1.31 Åpixel⁻¹ post-upgrade) on the red side, a compromise between broad spectral coverage and moderate spectral resolution. The plate scale was 043/pixel for both the blue and the new red detectors, and 078/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 (PA), 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 two or three separate exposures (pre- and 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.

Table 2. Observing Parameters

Object	Slit PA		Monitoring Period	Nobs	S/N	Nphot
	(deg)	(s)	(UT)
Zw 535−012	100	1800	2016/06/27−2017/03/04	41	57	10
I Zw 1	45	1600	2016/07/03−2017/02/01	34	103	6
Mrk 1048	0	1800	2016/08/08−2017/02/16	27	88	5
Ark 120	0	1800	2016/08/15−2017/03/26	34	124	5
Mrk 376	70	1800	2016/05/01−2017/05/01	31	84	6
Mrk 9	90	1800	2016/05/01−2017/05/01	33	78	5
Mrk 704	130	1800	2016/10/20−2017/05/01	23	106	3
MCG +04−22−042	145	1800	2016/05/01−2017/05/01	34	54	7
Mrk 110	44	1200	2016/05/01−2017/05/01	41	74	6
RBS 1303	20	1200	2016/05/01−2017/05/01	22	67	5
Mrk 684	60	1800	2016/05/01−2017/05/01	43	98	12
Mrk 841	50	800	2016/05/01−2017/05/01	45	77	11
Mrk 1392	50	1800	2016/05/01−2017/05/01	39	55	10
SBS 1518+593	90	1800	2016/05/01−2017/05/01	48	44	8
3C 382	70	1800	2016/05/01−2016/12/03	50	81	12
NPM1G+27.0587	60	1800	2016/05/01−2016/12/03	38	55	7
RXJ 2044.0+2833	60	1800	2016/05/01−2016/12/31	46	58	9
PG 2209+184	60	1800	2016/05/01−2016/12/31	40	32	9
PG 2214+139	60	1800	2016/05/18−2016/12/31	43	81	10
RBS 1917	30	1800	2016/06/01−2016/12/31	32	39	9
Mrk 315	60	1800	2016/06/07−2016/12/31	35	59	9

Note. 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.

Download table as:  ASCIITypeset image

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 six or fewer AGNs were observed when typically 12–15 objects would have been observed on a night with good conditions. We were confident that at least 14 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).

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 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 +04−22−042, and NPM1G+27.0587) and during particular nights where the seeing was exceptionally poor.

Most of the AGNs 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 AGNs 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 eight sources (Ark 120, RBS 1303, Mrk 9, RXJ 2044.0+2833, Zw 535−012, SBS 1518+593, Mrk 684, and Mrk 110), we further divided the [O iii] light curve by the median [O iii] flux, and subsequently normalized the Hβ light curve with the residual [O iii] light curve. This normalization helped to mitigate the residual scatter remaining in the Hβ light curves. We suspect that this flux-calibration anomaly due to the dichroic was occurring at a lower level throughout much of the campaign, which caused scatter larger than what might be typically achievable in other RM data sets.

For red-side data taken prior to 2016 September 20, each observation was flattened using a dome-flat exposure taken at the same telescope position. For red-side data taken between 2016 September 20 and 2017 January 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.

To ensure that all of 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 three to 12, with a median of eight 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 of 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 χ² 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 AGNs 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.

Zoom In Zoom Out Reset image size

Relative variability across the spectrum can be visualized using the rms spectrum. The Peterson et al. (2004) procedure of taking the standard deviation of flux values at each wavelength element over all of 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 broad-line 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 ii. The absence of 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.

Zoom In Zoom Out Reset image size

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 emission-line 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 AGNs and host-galaxy components. Residual errors from continuum subtraction can be particularly 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. 2013, 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) for the 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. 2021). 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 AGNs 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 spectral 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.

Zoom In Zoom Out Reset image size

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.

Table 3. Telescope and CCD Properties for Imaging Observations

Telescope	Dmirror	Detector	Field of View	Pixel Scale	Gain	Read Noise	Binning
	(m)			(″/pix)	(e−/ADU)	(e−)
FLWO-1.2 m	1.2	Fairchild CCD 486	231 × 231	0.336	4.45	7.18	2 × 2
KAIT	0.76	Apogee AP7	68 × 68	0.80	4.5	12.0	1 × 1
LCOGT-2 m	2.0	Merope	47 × 47	0.467	1.0	9.0	2 × 2
LCOGT-1 m	1.0	Fairchild Imaging	265 × 265	0.389	1.0	13.5	1 × 1
Liverpool	2.0	e2V CCD 231	10′×10′	0.304	1.62	8.0	2 × 2
MLO-1 m	1.0	Fairchild 446	133 × 133	0.718	2.13	3.88	1 × 1
Nickel	1.0	Loral	63 × 63	0.368	1.7	8.3	2 × 2
SPM-1.5 m	1.5	Fairchild 3041	54 × 54	0.32	4.2	14.0	2 × 2
WMO-0.9 m	0.9	Finger Lakes PL-09000	252 × 252	0.61	1.37	12.0	1 × 1

Download table as:  ASCIITypeset image

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 photon-counting 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 variations across the field of view. To account for these additional error sources, we measured the excess variance 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 , where σ_phot is the original photometric uncertainty on a given data point, and 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 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 cross-checked 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 Section 5.

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 (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 computed at every time lag τ at 0.5 day intervals via

where () 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³ realizations, building up distributions of correlation measurements. The median and ±1σ widths of the cross-correlation centroid distributions were then adopted as the final lag measurements, and their associated uncertainties as shown in Figures 4–10.

In any RM campaign, it is expected that some AGNs 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 AGNs, 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 AGNs (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. 2022; Li et al. 2021). Our method will be presented in detail by Guo et al. (in preparation), and we provide a brief description here.

We first generated light curves according to the damped random walk (DRW) model with the Python software CARMA ⁴² (Kelly et al. 2009, 2014). The DRW model provides an adequate description of ultraviolet/optical variabilities in AGNs, 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³ 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 of the simulations, using the same lag search range employed for each AGN, and counted, out of 10³, the number of positive lags (τ > 0) with peak values higher than our observed . 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 , 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 insignificant 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 , 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 and small values. Considering these lag-assessment results, we conclude that 16 of the 21 AGNs in our sample falling in the region and 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 compared to those selected by other large surveys (e.g., minimum r_max = 0.45 by SDSS-RM; Grier et al. 2017). These and values, along with AGN luminosities and lag measurements τ_cen in both observed and rest frames, are reported in Table 7. (Only the and are reported for the subsample of sources with unreliable lag measurements).

Zoom In Zoom Out Reset image size

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

Object	λ Lλ (5100 Å)	Observed Frame		Rest Frame				BLR Kinematics
		τcen	τpeak	τcen	τpeak
	(1043 erg s−1)	(days)	(days)	(days)	(days)
Zw 535−012	4.9 ± 0.7					0.63	0.12	Infalling
Mrk 1048	9.5 ± 1.8					0.62	0.11	Infalling
Ark 120	9.2 ± 3.4					0.95	0.03	Ambiguous
Mrk 9	5.6 ± 0.4					0.79	0.10	Ambiguous
Mrk 704	6.2 ± 0.6					0.82	0.19	Infalling
MCG +04−22−042	1.6 ± 0.4					0.95	0.00	Symmetric
Mrk 110	7.2 ± 1.7					0.92	0.02	Symmetric
RBS 1303	2.4 ± 0.4					0.90	0.02	Outflowing
Mrk 841	6.7 ± 0.9					0.74	0.02	Infalling
Mrk 1392	1.6 ± 0.6					0.95	0.01	Symmetric
SBS 1518+593	10.6 ± 1.0					0.60	0.03	Infalling
3C 382	15.0 ± 3.1					0.86	0.12	Symmetric
NPM1G +27.0587	10.9 ± 1.0					0.74	0.16	Infalling
RXJ 2044.0+2833	5.3 ± 0.5					0.78	0.04	Infalling
PG 2209+184	2.3 ± 0.6					0.89	0.01	Ambiguous
RBS 1917	4.9 ± 0.5					0.76	0.02	Outflowing
I Zw 1	28.5 ± 4.0	⋯	⋯	⋯	⋯	0.18	0.89	⋯
Mrk 376	17.6 ± 1.8	⋯	⋯	⋯	⋯	0.26	0.77	⋯
Mrk 684	8.3 ± 0.9	⋯	⋯	⋯	⋯	0.46	0.12	⋯
PG 2214+139	21.2 ± 1.3	⋯	⋯	⋯	⋯	0.62	0.45	⋯
Mrk 315	3.2 ± 0.4	⋯	⋯	⋯	⋯	0.53	0.23	⋯

Download table as:  ASCIITypeset image

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), and Grier 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 8. Velocity-resolved Lags

Object	Velocity Bin	τcen
	(km s−1)	(days)
Mrk110	−1910 ± 2324
Mrk110	−539 ± 417
Mrk110	−151 ± 357
Mrk110	205 ± 357
Mrk110	563 ± 357
Mrk110	1040 ± 596
Mrk110	1874 ± 1073

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

Download table as:  DataTypeset image

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 outflow-dominated 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 (L. Villafaña et al. 2021, in preparation).

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 and RBS 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 relates to other RM samples will be further examined in Section 6 and the Appendix, respectively.

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. 2009a, 2013; Kilerci Eser et al. 2015; Martínez-Aldama et al. 2020), but recent studies of increasingly 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. 2015, 2016; 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 radius–luminosity 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

with a median slope α = 0.45 and a median normalization K = 1.46 (Fonseca Alvarez et al. 2020). 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₀ = 67.8 km s⁻¹ Mpc⁻¹) provided by Wright (2006). All of our AGNs fall within 0.5 dex of the radius–luminosity relation as defined by the literature points within the range of λ L_λ (5100 Å) = 10^43.5–44.4 erg s⁻¹. The best-fit line in the form of Equation (5) for all of the data points combining our current work with those from the literature has a normalization K = 1.35 and α = 0.38.

Zoom In Zoom Out Reset image size

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 AGNs, the majority of Hβ velocity-delay maps appear roughly symmetric with rotationally dominated kinematics (e.g., Bentz et al. 2009b, 2010b; Denney et al. 2010; Barth et al. 2011; Horne et al. 2021). The AGNs with velocity-resolved reverberation detected to-date are mostly within the low-luminosity range [λ L_λ (5100 Å) = 10^42–43.5 erg s⁻¹]. Our results add significantly to the existing number of sources having velocity-resolved measurements with λ L_λ (5100 Å) ≈ 10⁴⁴ erg s⁻¹. 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 time-dependent 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 moderately high-luminosity sample appears consistent with that seen among AGNs with super-Eddington accretion rates (e.g., Du et al. 2016b).

Zoom In Zoom Out Reset image size

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 photon-counting 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 FWHM intensity of the line, and the line dispersion σ_line. The FWHM of a single-peaked line is defined as the difference between the wavelengths blueward and redward of the peak that correspond to half of its height. A double-peaked line requires additional attention to defining for the blue and red peaks, respectively, and thus their corresponding wavelengths at half the maximum height, but otherwise the procedure in determining the difference 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⁻¹ and dispersion σ_inst ≈162 km s⁻¹ for our small corrections. The resulting FWHM and σ_line measured from both the mean and the e-rms spectra are presented in Table 9.

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 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.

Table 10. Virial Products and Derived BH Masses

Object
	(107 M⊙)	(108 M⊙)
Zw 535−012
Mrk 1048
Ark 120
Mrk 9
Mrk 704
MCG +04−22−042
Mrk 110
RBS 1303
Mrk 841
Mrk 1392
SBS 1518+593
3C 382
NPM1G +27.0587
RXJ 2044.0+2833
PG 2209+184
RBS 1917

Note.^†Assuming f = 4.47 (Woo et al. 2015).

Download table as:  ASCIITypeset image

Ark 120 was monitored by Peterson et al. (1998) and Peterson et al. (2004) who had measured a rest-frame Hβ lag of 47.1 and days and σ_line of 1959 ± 109 and 1884 ± 48 km s⁻¹ from two different epochs, respectively. Adopting the same 〈f〉 of 4.47 from Woo et al. (2015) as in this study, the black hole masses calculated are 10⁸ M_⊙ and 10⁸ M_⊙, respectively.

In comparison, we measured a shorter lag of days in Hβ, resulting in a smaller black hole mass of 0.73 M_⊙. In a concurrent campaign, the Monitoring AGNs with Hβ Asymmetry (MAHA) group determined an Hβ lag of days and consequently a black hole mass of (using the same virial factor) in Ark 120 (Du et al. 2018a). Our measurement () is very similar to the MAHA result (), and roughly agrees with the previous masses () from over a decade ago within ±0.2 dex, indicating that the state of the AGN and the intrinsic BLR properties have stayed generally consistent in the past two decades.

Mrk 704 was among five bright Seyfert 1 galaxies that were targeted by De Rosa et al. (2018) in an RM campaign in early 2012. De Rosa et al. 2018 measured an Hβ lag of days and σ_line of km s⁻¹, resulting in a black hole mass of M_⊙, or . Four years later, we determined a longer lag of days and a larger σ_line of 2009 ± 102 km s⁻¹. Our black hole mass measurement of 1.03 M_⊙, or , is more than twice as large as the result from De Rosa et al. (2018), but our measurement’s uncertainty is substantial such that there is statistical agreement within 1.2σ.

Like Ark 120, Mrk 110 was another RM target monitored by Peterson et al. (1998) and Peterson et al. (2004) within the time period and had a rest-frame Hβ lag of 20.4, , and days with σ_line of 1115 ± 103, 1196 ± 141, and 755 ± 29 km s⁻¹ from three different epochs, respectively. Adopting the same 〈f〉 of 4.47 from Woo et al. (2015) as in this study, the black hole mass calculated correspond to 2.2 10⁷ M_⊙, 3.0 10⁷ M_⊙, and 1.7 10⁷ M_⊙, or a combined value of 2.3 10⁷ M_⊙ ().

With a measured Hβ lag of days and a σ_line of 1314 ± 69 km s⁻¹, we obtained a larger black hole mass of 3.5 M_⊙ (), consistent with the previous estimate within statistical uncertainties (1σ) and ± 0.2 dex.

SBS 1518+593 was one target from the MAHA campaign carried out by Du et al. (2018a) during a similar time frame as our LAMP study. Du et al. (2018a) established an Hβ lag of days and σ_line of 1038 ± 20 km s⁻¹. The inferred black hole mass adjusting for 〈f〉 is M_⊙, or . Measuring a similar lag of days and larger σ_line of 2249 ± 255 km s⁻¹, our black hole mass for this AGN was 5.5 M_⊙, or . Our measurements agree statistically within 1.3σ given the large uncertainties, with a difference roughly within ±0.2 dex. However, our velocity-resolved lag spectra disagree qualitatively. Du et al. (2018a) found the BLR kinematics to be in virialized motion while we see infalling behavior. This discrepancy may be due to differences in our respective campaigns’ monitoring durations, where they observed SBS 1518+593 for approximately 180 days while our coverage spanned 360 days. This discrepancy may also be due to the fact that our Hβ light curve exhibits much noise and scatter, which may in turn affect our velocity-resolved measurement.

The black hole mass in Zw 535-012 was previously determined to be 2.6 × 10⁷ M_⊙ (Wang et al. 2009), or . It was computed using the scaling law from Greene & Ho (2005), assuming (Hβ) = 42.1 erg s⁻¹ and FWHM(Hβ) = 2555 km s⁻¹ from a single-epoch spectroscopic observation in 2006. In our campaign, we obtained a black hole mass of 3.8 M_⊙, or , which is 1.4 times larger than the previous estimate but consistent within statistical uncertainties (0.6σ) and ±0.2 dex.