Monitoring AGNs with Hβ Asymmetry. IV. First Reverberation Mapping Results of 14 Active Galactic Nuclei

We report first-time reverberation-mapping results for 14 active galactic nuclei (AGNs) from the ongoing Monitoring AGNs with Hβ Asymmetry campaign (MAHA). These results utilize optical spectra obtained with the Long Slit Spectrograph on the Wyoming Infrared 2.3 m Telescope between 2017 November and 2023 May. MAHA combines long-duration monitoring with high cadence. We report results from multiple observing seasons for nine of the 14 objects. These results include Hβ time lags, supermassive black hole masses, and velocity-resolved time lags. The velocity-resolved lags allow us to investigate the kinematics of the broad-line region.


Introduction
Reverberation mapping (RM) is a spectroscopic monitoring technique that measures the masses of the supermassive black holes (SMBHs) fueling active galactic nuclei (AGNs) (e.g., Blandford & McKee 1982;Peterson 1993;Cackett et al. 2021).RM takes advantage of two defining observational signatures of AGNs: (1) variability on short time scales (days to months) and (2) broad emission lines (BELs) in their UV and optical spectra.These BELs arise from gas in the broad-line region (BLR), photoionized by the ionizing continuum.The BLR flux echoes the ionizing continuum flux-density variations some time later due to the finite speed of light.RM uses spectroscopic monitoring to measure the delayed response of the BLR, τ BLR , which is converted to an average radius of the BLR, R BLR , by multiplying by the speed of light.Assuming virialized motion of the BLR gas, we can use τ BLR along with the BLR velocity, measured from the width of a BEL, to estimate the SMBH mass.
Through great observational effort over the past three decades, more than 100 AGNs' SMBH masses have been measured with RM (e.g., Kaspi et al. 2000;Peterson et al. 2004;Kaspi et al. 2007;Bentz et al. 2009;Denney et al. 2010;Xiao et al. 2011;Grier et al. 2012;Barth et al. 2015;Edelson et al. 2015;Shen et al. 2015;Fausnaugh et al. 2016;Grier et al. 2017;Du et al. 2018a;Lu et al. 2022;Pandey et al. 2022;Chen et al. 2023; Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.Malik et al. 2023).RM anchors the empirical relationship between the BLR radius and AGN luminosity (Kaspi et al. 2000;Bentz et al. 2013;Du & Wang 2019) that allows for singleepoch mass estimation.
In addition to determining SMBH masses, RM can also reveal information about the BLR geometry and kinematics.With sufficiently high-cadence and signal-to-noise (S/N) data, we can measure τ BLR in velocity bins across an emission line.This allows for broad conclusions as to the dominant motion of the BLR gas (e.g., disklike, outflow, or inflow; Horne et al. 2004).There are ∼45 objects with velocity-resolved time lags (Bentz et al. 2008(Bentz et al. , 2009;;Denney et al. 2009Denney et al. , 2010;;Valenti et al. 2015;Du et al. 2016;Lu et al. 2016;Du et al. 2018b;De Rosa et al. 2018;Brotherton et al. 2020;Hu et al. 2020;Bentz et al. 2021;Bao et al. 2022;Bentz et al. 2022;Li et al. 2022;U et al. 2022;Bentz et al. 2023).Advanced techniques such as two-dimensional velocity-delay maps (Horne et al. 2004) and dynamical modeling of the BLR (Li et al. 2013;Pancoast et al. 2014) can reveal even more information about the BLR and provide a check on the fidelity of RM masses.
Monitoring AGNs with Hβ Asymmetry (MAHA) is an ongoing long-term RM campaign of objects with asymmetric Hβ emission lines.MAHA began in 2016 and has observed ∼100 objects for one or more seasons at different priority levels.Our primary targets have high-cadence, long-duration light curves that allow for not only SMBH mass estimations, but also advanced data analysis techniques such as measuring the line response in multiple velocity bins or dynamical modeling of the BLR gas.For more information on the MAHA campaign, see the first three MAHA papers, hereafter referred to as MAHA I (Du et al. 2018b), MAHA II (Brotherton et al. 2020), and MAHA III (Bao et al. 2022) as well as the following papers utilizing MAHA data: Kara et al. (2021), Oknyansky et al. (2021), Chen et al. (2023), and Oknyansky et al. (2023).We provide here a short summary of published MAHA results: MAHA I presented time lags, corresponding black hole masses, and velocity-resolved time lags for 3C 120, Ark 120, Mrk 6, and SBS 1518+593; MAHA II presented time lags, corresponding black hole masses, and velocity-resolved time lags for Mrk 79, NGC 3227, and Mrk 841.MAHA II also presented velocity-delay maps for Mrk 79 and NGC 3227.MAHA III presented time lags, corresponding black hole masses, and velocity-resolved time lags for 15 Palomar-Green quasars.In this paper, we present 14 MAHA objects without prior RM results.There are good results for multiple seasons for nine of the objects, and we show all the seasons here.We present velocity-resolved lags for all objects and discuss implications for the kinematics of the individual BLRs.
Section 2 describes our observations and data reduction as well as the construction of our light curves.Section 3 describes our analysis including time-lag measurements, line-width measurements, black hole mass calculations, and velocityresolved time-lag measurements.Section 4 contains the discussion of the results for each object.Finally, Section 5 summarizes our results.

Observations and Data Reduction
In this paper, we present results from RM of 14 MAHA AGNs that have not previously had successful time lags published.For more information about the MAHA target list as a whole, including sample selection criterion, see MAHA I, although we note that our targets have evolved over time.In short, we prefer bright, northern objects with asymmetric Hβ profiles as measured by the dimensionless asymmetry parameter where λ c (3/4) and λ c (1/4) are the central wavelengths at threequarters and one-quarter of the line peak and Δλ(1/2) is the FWHM of the line (De Robertis 1985; Boroson & Green 1992;Brotherton 1996).Blue asymmetries are positive and red are negative.Usually targets have been selected based on published or publicly available spectra showing asymmetric Hβ, but the profile may have varied by the time we began observing it.Target information, including coordinates and redshifts from NED,23 are listed in Table 1.We next describe our observations and data reduction.

Spectrophotometry
We conducted our observations using the Wyoming Infrared Observatory's (WIRO) 2.3 m telescope with the Long Slit Spectrograph primarily in remote operation mode (Findlay et al. 2016).We used a 900 line mm −1 grating which provides a wavelength range of 4000-7000 A ̊and a dispersion of 1.49A ̊pixel −1 .We used a 5″ wide slit oriented north-south.We observed at least one spectrophotometric standard star each night for flux calibration (either BD+28d4211, Feige 34, G191B2B, or HZ 44).We observed CuAr lamps for wavelength calibration and a uniformly lit white square for flat fields.We reduced the data with IRAF v2.16 using standard techniques and a uniform aperture of ±6 8 and background windows of 7 6-15 1 on both sides of the objects' nuclear emission.
There are sometimes flux-density variations between spectra due to clouds and slit losses.In order to obtain accurate photometry, we performed additional flux calibrations using the [O III] λ5007 line.This [O III] line arises from the narrowline region of the AGN which varies slowly, on the order of tens of years or longer, making it appropriate for use in an observation campaign of months (Peterson et al. 2013).We used the calibration technique of van Groningen & Wanders (1992) with minor modifications.Briefly, our wide slit (5″) is larger than the average FWHM of the seeing of our observations (2″-3″), and variable seeing between exposures leads to changes in the line point-spread function (PSF).In addition, telescope tracking errors can lead the line PSF to deviate from a Gaussian.To account for these effects we convolved the [O III] λλ4959, 5007 doublet profiles with a double Gaussian that fits the broadest [O III] profile observed for each object.This smooths the spectra and minimizes the effect of variable seeing and tracking inaccuracies.We then scaled each exposure of an object to match its standard O [III] λ5007 flux.Standard O [III] λ5007 fluxes are obtained from spectra taken in photometric conditions.For more details, see discussion in MAHA I. Our standard [O III] λ5007 fluxes are listed in Table 2.
We produced the final flux-calibrated spectra by averaging the exposures (appropriately noise weighted) from a single night.See Figure 1 for examples of typical MAHA spectra.
In order to investigate the general Hβ profile and measure line asymmetries and widths, we calculated mean and rms spectra using respectively.We present mean and rms spectra for each object in Figure 2. Note that the [O III] lines have essentially disappeared in the rms spectra, indicating data uniformity and consistency in our flux calibration.

WIRO Light Curves
Before constructing our light curves, we corrected each epoch of our calibrated spectra for Galactic reddening.We used the Cardelli et al. (1989) extinction law implemented by the Python package dust_extinction and the Schlafly & Finkbeiner (2011) R V = 3.1 and A V values for each object taken from NED.
We measured the 5100 A ̊continuum flux densities from our calibrated spectra by averaging the flux density spanning the wavelength range of 5075-5125 A ̊in the rest frame.We note that we do not remove host-galaxy starlight from our continuum flux densities.
There are two commonly used ways to measure emissionline fluxes: direct integration (e.g., Peterson et al. 1998a;Kaspi et al. 2000;Bentz et al. 2009;Grier et al. 2012;Du et al. 2014) or spectral fitting (e.g., Barth et al. 2013;Hu et al. 2015).Direct integration involves fitting a linear continuum underneath the Hβ profile, subtracting that continuum from the line profile, then integrating across the line profile.Spectral fitting involves detailed fitting of the various components of the spectrum to isolate the BEL.We followed MAHA I in using direct integration.To fit a linear continuum underneath the Hβ profile, we linearly interpolated across two emission-line-free windows on either side of the line profile.We used the rms spectrum to select the emission-line-free windows, choosing minimally varying windows to the left and right of the Hβ emission line.We also used the rms spectrum to select a Hβ   In addition to narrow-line contamination, the broad Hβ emission can also contain flux contamination from Fe II emission.Only two of our 14 objects show moderate optical Fe II emission: Mrk 813 and NGC 7603.Fe II varies on similar time scales as Hβ (Hu et al. 2015), consequently it should not strongly affect the Hβ lag measurement.Fe II emission is also relatively weak in the Hβ region, though it is stronger in the region of the O [III] doublet, making Fe II contamination a  2).We tried changing the Hβ integration limits for NGC 7603 in order to remove the O [III], and as a result the Fe II, emission from the Hβ integration range.This results in a lag measurement that is consistent with the lag measured using the full Hβ integration range, to within the uncertainties (see Table 2 for the Hβ integration range).This indicates that in NGC 7603 the Fe II emission that is blended with the broad Hβ emission does not significantly affect the time lag.We therefore elected to keep our full Hβ integration range for NGC 7603 and ignore the potential Fe II contamination.
The flux density and flux uncertainties in the light curve contain two components: (1) Poisson noise, and (2) random uncertainty from tracking error and variable observing conditions in different exposures of an object during the same night.We estimated the random uncertainty using the scatter of the flux density between consecutive exposures on the same night between ∼4100 and 5700 A ̊. We added these two sources of uncertainty in quadrature.We emphasize here that our error budget does not include the systematic uncertainty from our narrow-line model assumptions used to subtract the narrow-line contribution from Hβ.We also do not include uncertainty contributed by the narrow-line subtraction itself.Therefore our total Hβ flux uncertainties are likely underestimated.We plotted the total uncertainty as error bars in our light curves.

Supplemental Continuum Light Curves from Zwicky Transient Facility Photometry
We followed MAHA I in supplementing our continuum light curves with survey photometry when available and of shows the rms spectra for each analyzed season in units of 10 −16 erg s −1 cm −2 A ̊−1 .Both the mean and rms spectra are shifted relative to each other in order to display all seasons on one plot.The y-axis ticks and tick labels on the mean and rms spectra plots correspond to the maximum flux-density value of each spectrum shown.The blue region is the Hβ integration range and the gray regions are the line-free continuum regions used to fit the linear continuum under Hβ.All wavelengths are in the rest frame.The colors in the right panels correspond to the different seasons as defined in the middle-right panel.(The complete figure set (14 images) is available.)sufficiently small uncertainties.Contributions from optical emission lines (other than Hα) contaminating broadband filters is generally weak, making broadband photometry usually acceptable for continuum light curves.The additional data points provide improved temporal coverage and cadence.
The Zwicky Transient Facility is an all-northern-sky survey using the Palomar 48 inch Schmidt Telescope at Palomar Observatory.The survey goes to median depths of g ∼ 20.8 and r ∼ 20.6 mag.The survey began in 2018 March and builds upon the predecessor survey, the Palomar Transient Factory (Masci et al. 2019;IRSA 2022).When available, we used Zwicky Transient Facility (ZTF) g-band photometry to supplement our 5100 A ̊light curves.We chose g-band over r-band photometry as the g band contains (rest-frame) 5100 A ̊for objects with z 0.11, making it the appropriate choice for our sample (Bellm et al. 2019).We utilized the ztfquery package (Rigault 2018) in fetching the ZTF light curves.
It is necessary to scale the ZTF photometry to account for differing host-galaxy flux-density contributions due to aperture size differences between ZTF and WIRO.The scaling is a simple linear shift following To find the optimal a and b parameters we used the Bayesian-statistics-based package PyCali;24 PyCali assumes the light curve can be described by a damped random walk model (Li et al. 2014).We then used the best-fit parameters to scale the ZTF flux densities.We note that PyCali takes into account the uncertainty on a and b in the final uncertainty of the scaled ZTF flux densities.Finally, we combined data points from the same night using a weighted average.Visual inspection of the light curves show that the intercalibrated continuum data agree well.
See Table 3 for light-curve statistics, including the flux density and flux variation amplitude and uncertainty defined by Rodríguez-Pascual et al. (1997) as respectively.See WIRO light-curve data in Table 4, and object light curves in Figure Set 2 with the first object light curve shown in Figure 2.

Time-series Analysis
We measured the delay between 5100 A ̊continuum variations and Hβ line response using the interpolated cross correlation function (ICCF; Gaskell & Sparke 1986;Gaskell & Peterson 1987;White & Peterson 1994).The ICCF measures the strength of the correlation between the two light curves at a range of lags.It does this twice, alternating which light curve is interpolated.The output cross correlation function (CCF) is the average of the two.We can determine the time delay from the ICCF in two ways: the lag corresponding to the the peak correlation coefficient, r, or the centroid of some range of correlation coefficients (generally > 0.8 r max where r max is the maximum correlation coefficient).If the BLR has an extended geometry, the CCF peak is biased toward the inner radius while the centroid is a less biased, luminosity-weighted measurement (Koratkar & Gaskell 1991;Perez et al. 1992).The centroid is also more consistent with the virialized motion of the BLR gas (White & Peterson 1994).
To investigate time-lag uncertainties we used the flux randomization/random subset sampling (FR/RSS) method (Peterson et al. 1998b).The FR/RSS method accounts for measurement and sampling uncertainties.We generated mock light curves by perturbing the measured flux density or flux values according to their uncertainties and randomly resampling our light curves with repetition.The FR accounts for measurement uncertainties and the RSS accounts for uncertainty due to inclusion/exclusion of any one point.We used the FR/RSS method in a Monte Carlo fashion, running it 1000 times to build up cross correlation centroid (peak) distribution (CCCD or CCPD).We take the median and 68% confidence interval of the CCCD as our lag and uncertainty.We report the median of the CCPD with uncertainties for completeness.
We chose the lag search range for the CCF for each object based on visual inspection of the CCFs and CCCDs, choosing the upper limit that allows the CCF to reach a minimum.Our search windows all begin at −10 days and increase up to 105 days.Welsh (1999) showed that the CCF measurements can be strongly affected when a light curve has a varying mean or variance across the time sampled.In other words, the presence of long-term secular trends in a light curve will bias the timelag measurement toward smaller lags (see also Perez et al. 1992).The process of fitting and subtracting a low-order polynomial to the light curves is called detrending and has been explored in various RM studies (e.g., Grier et al. 2008;Peterson et al. 2014).We note that detrending might be removing long-term variation that is part of the reverberation response.We therefore proceed with caution and only detrend one light curve (see the discussion in Section 4.1) where there is a clear slope in the continuum light curve while the Hβ light curve stays mostly flat.Our detrending process consists of fitting a line to the light curve using nonlinear least-squares minimization and then subtracting that line from the light curve.We report time lags and lag search ranges in Table 6.

Line-width Measurements
The Hβ line is broad and can be contaminated by the narrow [O III] requiring a narrow-line subtraction in the mean spectra to accurately measure the line width.It is not necessary to subtract narrow lines from the rms spectra, as the narrow lines do not vary on RM time scales (Peterson et al. 2013).We followed the procedure outlined in MAHA I, fitting [O III] λ5007 as a single Gaussian and using the best-fit parameters as a template to fit [O III] λ4959.We assumed that the narrow component of Hβ has the same profile as [O III] and used our template to fit.We fixed the separation between  5 for the narrow-line flux ratios for each season of each object.
We used two methods to measure line widths: the FWHM and the line dispersion, σ Hβ .As many of our objects have significantly asymmetric or double-peaked Hβ profiles, we followed the procedure in Peterson et al. (2004) to measure the FWHM of double-peaked profiles, by considering the blue peak and the red peak separately.The second moment of the line profile is the mean square dispersion of the line We take the square root of the second moment as σ Hβ .
To evaluate the uncertainty in our line widths we follow the procedure established in MAHA I of using the bootstrap method.We resampled the mean/rms spectrum for each object, using N randomly selected points with repetition, out of the N points in the spectrum, where N in both cases is the number of points in the spectrum.We then measured the FWHM and σ Hβ from the resampled spectrum.We did this 500 times and took the median and standard deviation of the resulting distribution as our line-width measurement and uncertainty.
Finally, we subtracted instrumental line broadening from our line widths.MAHA I estimates the instrumental line-spread function by comparing the [O III] widths to the intrinsic widths reported in Whittle (1992).Following MAHA I, we adopted 925 km s −1 as the average line-spread function and subtracted it in quadrature from all line widths.This correction introduces a small uncertainty into the line-width measurements.We report line widths in Table 6

Black Hole Masses
The mass of the supermassive black hole fueling the AGN can be estimated assuming virialized motion (Peterson & Wandel 1999;Wandel et al. 1999) in the BLR using where R BLR = cτ BLR and ΔV 2 can be measured from the width of the Hβ line and f BLR is a calibration factor that accounts for the unknown geometry of the BLR.It is likely that each object has an individual f BLR , but dynamical modeling is needed to  find the best-fit BLR geometry for a given object and calculate its individual f BLR (e.g., Li et al. 2013;Pancoast et al. 2014;Li et al. 2018).As none of our sample targets have individual f BLR values in the literature, we use an average f BLR that is calibrated to bring active galaxy masses in line with the quiescent masses according to the M • -σ * or the M • -M * quiescent galaxy relationships (e.g., Onken et al. 2004;Woo et al. 2010;Graham et al. 2011;Park et al. 2012;Grier et al. 2013;Ho & Kim 2014;Woo et al. 2015;Batiste et al. 2017).This average f BLR is appropriate for a population of AGNs, but introduces an uncertainty of a factor of a few in individual objects.Following previous MAHA papers, we adopted the empirically calibrated average value from Woo et al. (2015) of f BLR 1.12 (4.47) when using FWHM (σ Hβ ) from the rms spectra.We do not include the uncertainty on f BLR in our SMBH mass-uncertainty calculations.Some care must be taken in the choice of lag and line-width measurement.As discussed in Section 3.1, we chose the median of the CCCD as our lag measurement.Regarding the line measurement, there are two choices to be made: mean or rms spectrum and FWHM or σ Hβ .The advantage of the rms spectrum is that it only displays the variable part of the spectrum, meaning that narrow lines and other nonvarying components drop out and we do not have to worry about them contaminating the broad-line-width measurement.This removes the need to do spectral fitting which can be ambiguous.Also, in accordance with the assumption of virialized motion in the BLR, we measured the lag of the BLR gas that is responding to continuum flux variations and the width of the rms spectrum measures the velocity of this same gas.The rms spectrum can be noisy, however, and it can sometimes be difficult to measure the line width.It also can be contaminated by other variable broad lines: He II λ4686 and Fe II in the case of Hβ.Also, the rms requires multiple epochs of spectra to construct, rendering it impossible for use in singleepoch mass estimations.The mean spectrum is generally clean and easily replaced with a single-epoch spectrum.It does require narrow-line fitting and subtraction which can be difficult if the red wing of Hβ blends with [O III] λλ4959 5007.For these reasons we measured line widths from both the mean and rms spectra, though we prefer the rms spectra because it is intuitive to measure velocities from the variable part of the line.
Regarding FWHM versus σ Hβ , it is important to note they are not necessarily interchangeable and have their own strengths and weaknesses (Peterson et al. 2004;Wang et al. 2020b).In practice, σ Hβ better follows the virial assumption, though it is sensitive to the line wings which can be blended with other features.FWHM on the other hand is more easily and precisely measured (though that is not true for the rms spectra which can be noisy; Peterson et al. 2004).See Dalla Note.The units for the continuum flux densities and uncertainties and Hβ fluxes and uncertainties are 10 −15 erg s −1 cm −2 A ̊−1 and 10 −13 erg s −1 cm −2 , respectively.
(This table is available in its entirety in machine-readable form.)Though we report both σ Hβ and FWHM line widths from both the mean and rms spectra, we only report SMBH masses using σ Hβ and FWHM measured from the rms spectrum.For completeness and to aid comparison with single-epoch black hole mass estimations, we also report simple virial products (assuming f Hβ = 1) measured using the FWHM from the mean spectrum.We report our masses and virial products in Table 7.We indicate our preferred SMBH mass values for each object with a "✓" in the last column of Table 7.We prefer SMBH masses that are calculated using σ Hβ in the rms spectra and have the smallest measurement uncertainties.
All our reported masses are individually relatively small and collectively span a small range.This is because we are limited to nearby, typically low-luminosity objects observable at WIRO and with time lags short enough to resolve over an observing campaign of several months.In addition, this paper focuses on MAHA objects that provide a result in a single season of observations.Future MAHA papers will report results from objects that require analyzing several observing seasons together to measure a time lag.

Velocity-resolved Time Lags
It is possible with high-S/N, high-cadence RM data to determine average time lags as a function of velocity across the line profile.Traditional RM recovers the average time lag for the observed season, but does not reveal information about the kinematics or geometry of the BLR.In the absence of intensive forward modeling, or potentially difficult to interpret velocitydelay maps, a reasonable method to estimate basic BLR kinematics is to determine the time lag as a function of velocity.MAHA is a high-cadence, long-term RM campaign and while dynamical modeling techniques will be used on MAHA data in the future, it is beyond the scope of the present effort.Instead, we measured velocity-resolved time lags and are able to make broad statements about the kinematics of the BLR for the objects presented here.Generally, the signatures of a disklike, outflow, or inflow BLR are as follows: disklike structures are indicated by shorter lags in the high-velocity gas and longer lags in the low-velocity gas; outflow shows "blue leads red" which means longer lags in the redshifted gas and shorter lags in the blueshifted gas; inflow shows "red leads blue" which means shorter lags in the redshifted gas and longer lags in the blueshifted gas.We note recent work suggesting that sometimes the outflow signature in velocity-resolved time lags can unexpectedly be "red leads blue" (Mangham et al. 2017).The BLR kinematics elude an easy description and we advise readers to take our interpretations based on velocity-resolved time lags as possibilities rather than certainties in most cases.
We used two different binning methods to calculate time lags as a function of velocity: equal velocity bins and equal rms flux-density bins.The advantage of equal rms flux-density bins is that they allot similar noise to each bin.We present both binning methods, as they should roughly agree when robust results are obtained.We divided each object's Hβ profile into as many bins as fit across the line, while constraining the bins to remain larger than our effective spectral resolution of 925 km s −1 .This ranges from 6 to 15 bins.We remeasured the binned flux densities using the integration method for each 1ES 0206+522 1 5.6  1.6 2.0 24.5  8.2 9.4 8.9  2.5 3.0 2 11.9  2.1 1.4 10.0  1.6 1.1 2.9  2.5 1.9 Mrk 813 4 22.9  8.7

✓
Note.The virial product uses line-width measurements from the mean profiles.The mass measurements use line widths measured from the rms profiles.The check marks in column (6) denote our preferred mass measurement for an individual object.We prefer the mass measurement using the s line 2 measured from the rms profile (column (5)).
(This table is available in machine-readable form.)epoch and then cross correlated each bin with the continuum, building up a CCCD for each bin.We used non-narrow-linesubtracted spectra for the velocity-resolved analysis.The velocity-resolved time lags are plotted in Figure 3. Table 8 lists the number of bins, kinematics of the BLR, and line asymmetry, with further discussion on individual objects in Section 4.

Discussion
Below we discuss results for individual objects.

1ES 0206+522
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.1.We selected this object for MAHA based on its asymmetric Hβ profile observed in 2012 (WIRO, private communication).Although this object showed a significantly less asymmetric profile at the beginning of MAHA, we included it in the campaign as it is bright, observable nearly year-round at WIRO, and historically showed asymmetric Hβ.We observed this object for five seasons beginning in 2018.In the most recent season (season 5), 1ES 0206+522 started to show stronger blue asymmetry.We measured a time lag in three of those seasons, albeit changing lags across all three seasons, with season 1 showing the shortest and season 2 the longest lag.Season 1 has the largest mass measurement (8.9 × 10 7 M e ) and season 2 has the smallest (4.1 × 10 7 M e , which is our preferred mass).We detrended the season 1 continuum light curve, noting that the continuum light curve shows a long-term decreasing trend while the Hβ light curve remains mostly flat until the last few epochs of the season.The time lag, CCF, and CCCD for season 1 are measured from the detrended light curve.The velocityresolved lags indicate changing kinematics in different seasons.Seasons 1 and 5 suggest a disklike structure.Season 2 has a nearly constant lag across all velocities (see Figure 3).We excluded seasons 3 and 4 as we did not observe a clear reverberation.

Mrk 1040
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.2.Mrk 1040 is a bright Seyfert galaxy without particularly strong Hβ asymmetry in our first season, and was initially selected to fill a gap in the observing schedule near the Galactic plane during long winter nights.We dropped Mrk 1040 after one season to make room for targets with more asymmetric lines, but later resumed to obtain a second epoch of reverberation mapping to look for any changes in the Hβ profile and any corresponding changes in the BLR kinematics.Seasons 1 and 3 provide time lags, while season 2 is inconclusive.The CCF for season 3 is much broader than for season 1, showing a large "shoulder" for lags longer than the measured lag.This is likely because the continuum light curve shows two inflection points, and the Hβ light curve could be matched up with either inflection point.The mass measured in season 3 is larger than that measured in season 1 (2.1 × 10 7 and 7.0 × 10 6 M e , respectively).The season 1 mass is our preferred mass.Aside from the first bin, the velocity-resolved lags for season 1 show a general disklike pattern.The velocity-resolved lags for season 3 indicate outflow (see Figure 3).The asymmetry is red in both seasons, though it is stronger in season 3.

Mrk 618
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.3.Mrk 618 is a bright Seyfert galaxy with asymmetric Hβ and is notably a potential GRAVITY target given its southern decl.(Wang et al. 2020a).Mrk 618 was observed as part of a 2012 De Rosa et al. (2018) RM campaign, but they did not obtain a robust measurement of an Hβ time lag.We have been observing this object since 2019 and have four seasons of data.We measure a time lag from all four seasons.Different time lags are measured for the four seasons, ranging from 9.2 to 30.9 days.Season 2 displays the strongest variability and clearest reverberation.The CCF for season 3 is comparatively broad and shows two peaks, making the season 3 lag measurement the most uncertain.The black hole masses from seasons 1 and 2 and 1 and 4 agree to within the measurement uncertainty.The season 4 mass is our preferred mass.The velocity-resolved lags for season 2 looks generally disklike, while season 4's kinematics indicate outflow.The velocity-resolved lags for seasons 1 and 3 are indeterminate (see Figure 3).All four seasons show red asymmetry with season 1 showing the strongest asymmetry.

MCG -02-14-009
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.4.This target is another bright southern Seyfert galaxy that may be a good target for GRAVITY (Wang et al. 2020a).We have been observing MCG -02-14-009 since 2019 and have five seasons of data.Seasons 1 and 2 are nonresults due to a low number of epochs.We measured similar time lags for seasons 3 and 4 (12.1 and 11.5 days, respectively) and a longer lag for season 5 (17.4 days).The uncertainties on these lags are large and all three agree to within the measurement uncertainties.The season 5 CCF is double peaked with a broad CCCD.As a result, the season 5 lag is the most uncertain.The three mass measurements agree to within the margin of error.The season 3 mass is our preferred mass.The velocity-resolved time lags are indeterminate for all three seasons (see Figure 3).All three seasons show red asymmetry, with season 5 having the strongest red asymmetry.

IRAS 05589+2828
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.5.We initially targeted this object as its Hβ was described as having "a red wing" (Winter et al. 2010).This object is part of a close dual AGN with companion 2MASX J06021107+2828382 at a separation of 8 kpc (∼12″; Koss et al. 2012;Benítez et al. 2022).We have five seasons of data on IRAS 05589+2828, beginning in 2018.This object shows clear variation in each season and also shows a significant flux density and flux drop in season 3 with a sharp flux density and flux rise in season 4. We are able to measure time lags for each season.The lags are similar in seasons 1-3, about 20 days, and increase in seasons 4 and 5 to about 40 and 33 days, respectively.The mass measurements from seasons 1-5 agree within the uncertainties.The season 3 mass is our preferred mass.The velocity-resolved results for this object show a disklike structure for season 1, and evidence of inflow for season 2. The velocity-resolved lags for Seasons 3-5 are indeterminate (see Figure 3).All seasons show red asymmetry with season 3 showing the strongest red asymmetry.

Mrk 715
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.6.We described this object, also known as SDSS J100447.61+144645.6, as having "a doublepeaked Hβ line profile and a long tail to the red" (Du et al. 2018b).We have six seasons of data on Mrk 715, but were only able to measure time lags in seasons 4 and 5. Seasons 1-3 and 6 did not have favorable variability and our temporal coverage was not extensive.The variable component of Hβ in the rms spectra of this object is weak and appears dominated by a variable continuum.We followed MAHA III in fitting an AGN power-law continuum to the spectra using a nonlinear leastsquares minimization fit and subtracting the AGN power-law Correspondingly, the season 4 mass is more than twice that of the season 5 mass given the similar line widths.The season 5 mass is our preferred mass as it has smaller measurement uncertainty.The velocity-resolved results for season 4 show outflow, while season 5 shows a mostly disklike structure, though the rms profiles are noisy and we caution against giving these velocity-resolved measurements too much weight (see Figure 3).Both seasons show red asymmetry with season 5 having the stronger asymmetry.

SBS 1136+594
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.7.This object was selected based on the asymmetric profile published by Zamfir et al. (2010), which displays bumps on both the blue and red wings, with excess emission on the red side.We have two seasons of data, both showing a strong inflection point in the continuum with a clear, delayed response in Hβ.In season 2 we measure a lag more than twice that in season 1, though they agree to within the uncertainties.The CCCD for season 2 is far broader than the CCCD for season 1, leading to relatively larger uncertainties in the season 2 lag and mass measurements.Season 1 is our preferred mass.The velocity-resolved results are indeterminate in season 1 and indicate outflow in season 2 (see Figure 3).Both seasons show Hβ having red asymmetry with season 1 having the stronger asymmetry.

VIII Zw 233
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.8.We described this object as having "a redshifted peak" (Du et al. 2018b).We have been observing this object for five seasons, since 2016.Only seasons 2 and 3 produced a measurable time lag.Season 2 shows an apparent flare in the continuum light curve and a strong response a short time later from Hβ.The variable component of Hβ in the rms spectra of this object is weak and appears dominated by a variable continuum.We follow MAHA III in fitting an AGN power-law continuum to the spectra using a nonlinear leastsquares minimization fit and subtracting the AGN power-law continuum from each epoch of spectra before constructing the rms spectra for both seasons 2 and 3.After calculating this continuum-subtracted rms spectrum, it is clear that season 2 has much stronger variability in Hβ than season 3. The season 2 mass is our preferred mass.The velocity-resolved lags show a disklike structure in both seasons (see Figure 3).Both seasons show blue asymmetry.

Mrk 813
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.9.This object is a bright Seyfert galaxy that shows a bump and excess emission on the red side of the Hβ profile in SDSS spectra, although the line appears more symmetric in our most recent spectra.We have four seasons of data on this object, but not every year as the priority was dropped after poor initial results.Reprioritizing Mrk 813 in 2022 resulted in finally measuring a time lag, although the situation is not totally straightforward.The CCF for season 4 is double peaked which leads to a broad CCCD.The two peaks are separated by approximately 40 days.The two peaks are not strongly different than a broad, singly peaked CCF, and so we trust our CCCD which indicates a lag in the center of the two peaks.This season shows red asymmetry.The velocityresolved results for this object show inflow (see Figure 3).4.10.SDSS J145307.92+255433.0.0 See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.10.This object has a strange Hβ profile in a Sloan Digital Sky Survey spectrum and also appeared double peaked and asymmetric in our WIRO data.We only have one season of observations for this object and ZTF photometry helps significantly in covering gaps in the continuum light curve.The variation in our single season is strong and echoed by the Hβ.This season shows red asymmetry.The velocity-resolved results for this object show a likely disklike structure (see Figure 3).

SDSS J152139.66+033729.2
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.11.We previously described this object as having "an Hβ line with a red asymmetry" (Du et al. 2018b).We have four seasons of this object, extending back to 2017, but only observed an inflection point in both the continuum and Hβ in season 3, allowing us to measure a time lag.The lag in season 3 is clear, and the velocity-resolved results indicate a disklike structure (see Figure 3).This season shows red asymmetry.

2MASX J21090996-0940147
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.12.This target is another bright southern object that may be a good target for GRAVITY (Wang et al. 2020a).We have four seasons of data on this object, beginning in 2019, but only measure a lag in three seasons (1 and 3-4).Our measured lag is shortest in season 1 and longest in season 4 (ranging from about 6 to 30 days), with masses following the same trend (ranging from 5.0 × 10 6 to 2.2 × 10 7 M e ).Though the lags change by quite a bit, there is not a significant luminosity change between seasons.Season 1 is our preferred mass.Seasons 1 and 4 show blue asymmetry while season 3 shows red asymmetry.The velocity-resolved time lags for season 3 are suggestive of a disk and are indeterminate for seasons 1 and 4 (see Figure 3).

PG 2304+042
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.13.The only Palomar-Green quasar in this paper, the Hβ profile is asymmetric with a bump on the red wing.We have three seasons of data on this object, starting in 2020.All three seasons display sufficient variability in order to measure a time delay and the delays agree for each season.The measured masses also agree well.The kinematics from the velocity-resolved results suggest a disk in season 1, are flat in season 2, and suggest outflow in season 3 (see Figure 3).All three seasons show red asymmetry with season 1 having the strongest asymmetry.

NGC 7603
See the light curve, CCFs, CCCDs, mean, and rms spectra for this object in Figure 2.14.This object, also known as Mrk 530, has a history of variable and asymmetric Hβ profiles (Goodrich 1989;Kollatschny et al. 2000).This object has moderate Fe II emission in its optical spectrum.See Section 2.2 for an explanation of how we addressed the possible Fe II emission contaminant when constructing our Hβ light curves.This object has one season of data.The velocity-resolved results indicate a disklike structure (see Figure 3).This season shows red asymmetry.

Conclusions
In this fourth paper in the MAHA series we report RM results, including velocity-resolved time lags, for 14 objects without prior RM results.We present multiple seasons of data for nine of the 14 objects, allowing us to investigate time lag and BLR kinematic evolution over a period of several years.We find evidence for changing time lags without corresponding luminosity changes as well as changing line asymmetries and BLR kinematics as revealed by velocity-resolved lags.We observe predominantly red asymmetries.We observe mostly disklike BLR kinematics, with 12 out of 33 observed seasons displaying disklike kinematics.Five out of the 33 seasons show outflow, two show inflow, and two show the same time lag across the line profile.Twelve out of the 33 have kinematics that are not determinable from their velocity-resolved time lags.
integration window that contains the line as shown in the rms spectrum.Table2 listsHβ integration windows for each object.The broad Hβ emission contains flux contamination from the narrow Hβ component and, in some objects, the [O III] λλ4959, 5007 doublet.We isolated the broad component of the Hβ emission line by subtracting any narrow-line component from our measured Hβ fluxes.We began by fitting the [O III] doublet and narrow Hβ in each individual spectrum (for the fitting details see Section 3.2).We then integrated the fits to obtain the narrow-line fluxes.Finally we subtracted from each epoch the narrow-line fluxes included in each object's Hβ profile as determined from their rms profiles.

Figure 1 .
Figure 1.Randomly selected typical WIRO spectra for each object analyzed in this paper.The spectra are reduced and flux calibrated as discussed in Section 2.2.Flux-density units are 10 −15 erg s −1 cm −2 A ̊−1 .Wavelengths are in the rest frame.The vertical dashed blue line indicates the center of the Hβ line.Each spectrum is labeled with the name of the object and civil date observed.

Figure 2 .
Figure 2. Time-series analysis of 1ES 0206+522.The top-left panel shows the WIRO continuum and Hβ light curves for each season in arbitrary flux-density or flux units, respectively, shifted relative to each other so as to both display on the plot.The grayed-out seasons do not give a reliable result.The middle-left panel shows combined WIRO and ZTF continuum light curves for each analyzed season in units of 10 −15 erg s −1 cm −2 A ̊−1 .The orange line in season 1 shows the detrending line fit to the continuum.The orange season 1 light curve shows the detrended light curve.The bottom-left panel shows Hβ light curves for each analyzed season in units of 10 −13 erg s −1 cm −2 .The top-right panel shows CCFs and CCCDs in the rest frame for each analyzed season (the season 1 CCF/CCCD is measured using the detrended continuum light curve).The middle-right panel shows the AGN power continuum-subtracted mean spectrum for each analyzed season in units of 10 −16 erg s −1 cm −2 A ̊−1 .The dashed line shows the [O III] doublet and Hβ narrow-line components and the solid line shows the Hβ broad component.The bottom-right panelshows the rms spectra for each analyzed season in units of 10 −16 erg s −1 cm −2 A ̊−1 .Both the mean and rms spectra are shifted relative to each other in order to display all seasons on one plot.The y-axis ticks and tick labels on the mean and rms spectra plots correspond to the maximum flux-density value of each spectrum shown.The blue region is the Hβ integration range and the gray regions are the line-free continuum regions used to fit the linear continuum under Hβ.All wavelengths are in the rest frame.The colors in the right panels correspond to the different seasons as defined in the middle-right panel.(The complete figure set (14 images) is available.) λ4959, and Hβ to their laboratory separations.The [O III] λ5007/[O III] λ4959 and [O III] λ5007/Hβ flux ratios measured from our line fits are close to typical AGN values, 2.98 and 10, respectively (Kewley et al. 2006; Stern & Laor 2013), indicating the fitting procedure outlined here works well.See Table ). Columns (3)-(5) are the duration, number of epochs, and average cadence of the spectroscopy in days.Columns (6) and (7) are the mean flux density in units of 10 −15 erg s −1 cm −2 A ̊−1 and the fluxdensity variance of the continuum light curve.Columns (8) and (9) are the mean flux in units of 10 −13 erg s −1 cm −2 and the flux variance of the Hβ light curve.The uncertainty on the mean flux densities and fluxes is the standard deviation of the light curves.

Figure 3 .
Figure3.Velocity-resolved lags.The top panel in each subplot shows the lags using bins of equal velocity width.The middle panel shows the lags using bins of equal rms flux density.In both the top and middle panels, the cross correlation coefficients for each bin are written at the top of the bin.The bottom panel shows the rms spectrum without continuum subtraction in units of 10 −16 erg s −1 cm −2 A ̊−1 .The bottom panel also shows the object name and season number.Grayed-out points are lag measurements with cross correlation coefficients less than 0.5.

Table 2
Measurement Windows and [O III] Standard Fluxes The second column contains the O [III] λ5007 standard fluxes used for flux calibration.The third and fourth columns contain the windows on the blue and red side of Hβ used for fitting a linear continuum underneath the line.The fifth column contains the window used for integrating across the Hβ profile.All wavelengths are in the rest frame.

Table 7
Virial Products and Black Hole Masses

Table 8
Velocity-resolved Results and Hβ Asymmetry Note.Column (1) lists the object name.Column (2) lists the season.Column (3) lists the number of velocity bins.Column (4) lists the BLR dynamics as indicated by the velocity-resolved time lags.Column (5) lists the asymmetry where negative numbers indicate red asymmetry and positive numbers indicate blue asymmetry.Column (6) lists the Hβ integration range when it is different than the range used for constructing the light curves in Table2.Wavelengths are in the rest frame.