Monitoring AGNs with H\beta\ Asymmetry. I. First Results: Velocity-resolved Reverberation Mapping

We have started a long-term reverberation mapping project using the Wyoming Infrared Observatory 2.3 meter telescope titled"Monitoring AGNs with H\beta\ Asymmetry"(MAHA). The motivations of the project are to explore the geometry and kinematics of the gas responsible for complex H\beta\ emission-line profiles, ideally leading to an understanding of the structures and origins of the broad-line region (BLR). Furthermore, such a project provides the opportunity to search for evidence of close binary supermassive black holes. We describe MAHA and report initial results from our first campaign, from December 2016 to May 2017, highlighting velocity-resolved time lags for four AGNs with asymmetric H\beta\ lines. We find that 3C 120, Ark 120, and Mrk 6 display complex features different from the simple signatures expected for pure outflow, inflow, or a Keplerian disk. While three of the objects have been previously reverberation mapped, including velocity-resolved time lags in the cases of 3C 120 and Mrk 6, we report a time lag and corresponding black hole mass measurement for SBS 1518+593 for the first time. Furthermore, SBS 1518+593, the least asymmetric of the four, does show velocity-resolved time lags characteristic of a Keplerian disk or virialized motion more generally. Also, the velocity-resolved time lags of 3C 120 have significantly changed since previously observed, indicating an evolution of its BLR structure. Future analyses of the data for these objects and others in MAHA will explore the full diversity of H\beta\ lines and the physics of AGN BLRs.


INTRODUCTION
The most prominent features in the UV and optical spectra of luminous active galactic nuclei (AGNs) and quasars are the broad emission lines (BELs), with velocity widths ranging from ∼ 10 3 km s −1 to ∼ 10 4 km s −1 (e.g., Schmidt 1963;Osterbrock & Mathews 1986;Boroson & Green 1992;Sulentic et al. 2000;Shen et al. 2011, and references therein). The spectra generally show similarities from local Seyfert galaxies all the way to z 7 quasars (e.g., Francis & Koratkar 1995;Vanden Berk et al. 2001;Mortlock et al. 2011;Bañados et al. 2018). This suggests that the broad-line regions (BLRs) of AGNs have similarities from object to object over cosmological time. The similarities lead to recogni- * PI of the MAHA Project tion of a common formation mechanism of the BLRs in most AGNs and quasars, and also to the importance of investigating nearby AGNs for insight into high−z quasars. However, we so far possess insufficient understanding of the BLRs and their physics.
Reverberation mapping (RM) campaigns (e.g., Peterson et al. 1993Peterson et al. , 1998Peterson et al. , 2002Kaspi et al. 2000Kaspi et al. , 2007Bentz et al. 2008Bentz et al. , 2009Denney et al. 2009;Barth et al. 2011Barth et al. , 2013Barth et al. , 2015Rafter et al. 2011Rafter et al. , 2013Du et al. 2014Du et al. , 2015Du et al. , 2016bDu et al. , 2018Wang et al. 2014;Shen et al. 2016b;Fausnaugh et al. 2017;Grier et al. 2012Grier et al. , 2017b provide measurements of BEL lags (τ BLR ) with respect to the varying ionizing continuum as originally suggested by Bahcall et al. (1972) and Blandford & McKee (1982). The time lag multiplied by the speed of light (c) provides the size scale of the BLR. In individual objects, different lines echoing at different distances are consistent with the same virial black hole mass indicating Keplerian motion, where G is the gravitational constant, R BLR = cτ BLR is the emissivity-weighted radius of the BLR, ∆V is either full-width-half-maximum (FWHM) or velocity dispersion (σ line ) of mean spectra or root-meansquare (rms) spectra Wandel et al. 1999), and f BLR is the virial factor determined by the geometry and kinematics of the BLR (e.g., Peterson et al. 2004).
In past decades, RM established a tight correlation between time lags and luminosity based on a heterogeneous AGN sample containing mostly sub-Eddington sources (Kaspi et al. 2000;Bentz et al. 2013). More recently, however, new RM campaigns have discovered that lags are substantially shorter for super-Eddington AGNs than their sub-Eddington counterparts of the same luminosity (Du et al. , 2016b(Du et al. , 2018. Beyond just determining a single global time lag, high-quality reverberation mapping can provide time lags as a function of velocity. Such velocity-resolved time lags are diagnostic of the BLR structure itself, and these and other advanced analysis products such as velocity-delay maps and dynamical models commonly find that the BLR is comprised of a thick Keplerian disk (e.g., Bentz et al. 2008;Denney et al. 2010;Grier et al. 2013b;Pancoast et al. 2014;Du et al. 2016a;Grier et al. 2017a) 1 . This flattened disk of BLRs as the major ingredient could be common in AGNs and quasars. The BLR, however, appears significantly more complex in many objects. The BEL profiles (especially those of the Hβ line) of most AGNs are roughly symmetric (e.g., Gaussian or Lorentzianlike), but there are still a large number of objects showing asymmetric (redward, blueward or even double-peaked) profiles both for low and high-ionization lines (e.g., Boroson & Green 1992;Eracleous & Halpern 1994;Brotherton 1996;Peterson 1997). Similarly, Hβ usually shows a line peak close to systemic velocity, but extreme outliers exist (Eracleous et al. 2012;Shen et al. 2016a). Some net radial motion and/or opacity effects would seem required to explain BLRs displaying extreme profiles. The BELs with the extreme profiles may indicate special geometry and kinematics of their BLRs, such as super-fast or ultra-strong inflow or outflow, or even abnormal nuclear environments (e.g., absorption or dust). The asymmetries of Hβ profiles significantly correlate with the ratio of Fe II to Hβ (Boroson & Green 1992), and with [O III] luminosities (Wang et al. 2017a). High-precision RM has shown indications of infall and (more rarely) outflow 1 There are several other lines of evidence indicate that the typical BLR has the structure and dynamics of a flattened disk with emission-line gas following Keplerian motion (see Gaskell 2009 for a review of the BLR). For instance, Hβ velocity widths are systematically smaller in more face-on systems as determined using orientation-dependent radio properties (Wills & Browne 1986), and spectropolarimetry of type 1 AGNs also supports a flattened, disk-like geometry (Smith et al. 2005;Baldi et al. 2016;Savić et al. 2018  in Hβ-emitting gas, sometimes in the presence of a Keplerian disk component as well (e.g., Grier et al. 2017a). Some objects with excellent RM data sets have defied simple explanation, such as Mrk 6 (Doroshenko et al. 2012;Grier et al. 2013b). Mrk 6 is also noteworthy for its complex Hβ profile, which possesses a strongly blueshifted peak in addition to a lower velocity peak in most epochs, and red asymmetric tail. More detailed investigations of such AGNs are needed to provide a deeper understanding of the full diversity of the structure and kinematics of the BLR, and perhaps clues to its origin as well. The latter has been a controversial subject that is not yet resolved, although there are many proposals (Murray & Chiang 1997;Czerny & Hryniewicz 2011;Elvis 2017;Wang et al. 2017a;Baskin & Laor 2018).
Another intriguing hypothesis to explain some asymmetric and shifted profiles is the existence of close binary supermassive black holes (CB-SMBHs) in AGNs. CB-SMBHs have been predicted to be located in galactic centers due to galaxy mergers (Begelman et al. 1980), however, one single black hole is commonly assumed in explanations of RM data. This assumption may be challenged by the more unusual profiles of BELs implying the potential presence of CB-SMBHs (e.g., Gaskell 1996;Eracleous et al. 1997Eracleous et al. , 2012Boroson & Lauer 2009;Bon et al. 2012Bon et al. , 2016Li et al. 2016Li et al. , 2017Decarli et al. 2013;Shen et al. 2013;Wang et al. 2017b;Runnoe et al. 2017;Nguyen & Bogdanović 2016;Pflueger et al. 2018), which could be an indicator of CB-SMBHs (Popovic et al. 2000;Shen & Loeb 2010). However, there are several alternative explanations for the complicated profiles, such as elliptical BLR disks for asymmetric or double-peaked profiles (Eracleous et al. 1997), hot spots (Jovanović et al. 2010), partially covering dusty obscurers (Gaskell & Harrington 2018), and even spiral-arms (Storchi-Bergmann et al. 2017). Recoiling AGNs may be an additional possibility if the velocity is high enough to create an appreciable shift from a narrow-line region reflecting the redshift of the host galaxy (Volonteri & Madau 2008).
CB-SMBHs are expected to be sources of nano-Hz gravitational waves, likely to be discovered by Pulsar Timing Arrays (e.g., Sesana et al. 2008 for AGNs that are binary candidates in order to provide the optical identification of CB-SMBHs. For these reasons, we have undertaken the "Monitoring AGNs with Hβ Asymmetry" (MAHA) project. We describe our target selection, observations, and data reduction in Section 2. We report the initial results for four AGNs here, specifically: the mean and rms spectra in Section 2.3; the measurement of the light curves in Section 3; the lag measurements, the widths of the Hβ lines, black hole masses, and the velocity-resolved time lags in Section 4. We discuss the results for each individual object in Section 5. Finally, we summarize our findings in Section 6. Below we describe our MAHA target selection as well as our program of observations and data reduction that we started in December 2016. Due to a lack of significant variability, long time lags comparable to the length of an individual seasonal campaign, or other issues, many objects will require multiple campaigns to produce high-quality results of the type we seek. We expect to add or drop objects as MAHA progresses for a variety of reasons, so the sample should not be considered as absolute. The sample should, however, illustrate the type of objects we are monitoring and the diversity of their Hβ profiles. Our observational methods and data reduction apply generally to the entire sample, although we provide measurements and analyses in this first paper for only four objects for which we have obtained good quality velocity-resolved time lags.

Targets
The core MAHA sample includes AGNs with asymmetric or double-peaked Hβ emission lines, as well as objects reported as binary black hole candidates (which usually also have asymmetric lines). While the majority of the Hβ profiles of MAHA targets can easily be visually identified as asymmetric, it will be useful to parameterize asymmetry. There exist many ways to quantify asymmetry, each with its pros and cons. For ease of historical comparison to previous work (e.g., De Robertis 1985; Boroson & Green 1992;Brotherton 1996), we adopt the dimensionless asymmetry parameter: where λ c (3/4) and λ c (1/4) are the central wavelengths at the 3/4 and 1/4 of the peak height, and ∆λ(1/2) is the FWHM. Blue asymmetries are positive, red are negative. Note that the asymmetry is independent of any peak or centroid wavelength shift relative to systemic. Figure 1 illustrates what we mean when we say red and blue asymmetries and doublepeaked lines. We selected our targets from a variety of literature sources based on both asymmetry measurements and visual inspection of optical spectra (e.g., De Robertis 1985;Stirpe 1990;Boroson & Green 1992;Eracleous & Halpern 1993, 1994Marziani et al. 2003;Hu et al. 2008a,b;Eracleous et al. 2012). We selected some additional sources directly from the quasar sample of the Sloan Digital Sky Survey (SDSS) (Schneider et al. 2010), decomposing the Hβ lines by two Gaussians through the multi-component fitting procedure in Hu et al. (2008a,b), then selecting objects for which the two Gaussians have substantially different central velocities. This method recovered many objects already identified from the literature sources above.
Using the Wyoming Infrared Observatory (WIRO) 2.3 meter telescope and its longslit spectrograph requires additional selection criteia. WIRO has a latitude of 41 • N, and we do not include targets south of declination ∼ −10 • . In order to obtain sufficiently high signal-to-noise ratios (S/N) in exposure times of ∼1 hour or less, the magnitude (in V or SDSS r band) needs to be brighter than ∼17. We also require z < 0.38 to keep [O III] at less than 7000Å observed frame to avoid the inefficiencies of grating tilts and extra calibrations. Finally, because of our procedure for photometric calibration using narrow lines, [O III]λ5007 cannot be too weak. When everything else is equal, we prefer brighter targets at lower redshifts and more northern declinations. We do not automatically exclude AGNs from MAHA just because they have previous RM results. BLRs and their corresponding Hβ asymmetries may evolve over time periods of several years. Additionally, we aim to obtain high enough data quality to not only determine velocity-resolved time lags, but also to conduct more involved analyses such as creating velocity-delay maps and dynamical models (e.g., Horne 1994;Horne et al. 2004;Bentz et al. 2010;Grier et al. 2013b;Skielboe et al. 2015;Pancoast et al. 2011Pancoast et al. , 2012Pancoast et al. , 2014Grier et al. 2017a).
We provide additional information about the MAHA targets and their spectra in Appendix A, including the initial MAHA sample, its characteristics, and example WIRO spectra of the Hβ line region obtained during 2016-2018. Because of the need to use [O III] to calibrate fluxes in WIRO spectra, and the fact that objects with strong [O III] tend to be the ones with the strongest red asymmetric Hβ lines (Boroson & Green 1992), our sample is biased against objects with blue asymmetric Hβ lines. We plan to search for and add more AGNs with blue asymmetric Hβ lines as MAHA progresses.
As the first paper of the series, we provide here our RM measurements of 4 AGNs: 3C 120, Ark 120, Mrk 6, and SBS 1518+593 (see Table 1 for their coordinates and some general information). Their luminosities place them among Seyfert galaxies and broad line radio galaxies rather than the more energetic quasar category. Ark 120 and Mrk 6 show extreme red asymmetric profiles, while 3C 120 and SBS 1518+593 show milder red asymmetries. None are extreme super-Eddington accretors, although 3C 120 and SBS 1518+593 have dimensionless accretion ratesṀ on order unity , see also 4.3). The light curves of all these sources have shown unambiguous peaks or dips, and are sufficient for us to measure Hβ time lags.

Spectroscopy
We obtained spectroscopic data using the 2.3 m telescope at the WIRO and its Long Slit Spectrograph, observing remotely from the University of Wyoming campus (Findlay et al. 2016). We used the 900 line mm −1 grating, which pro-vides a dispersion of 1.49Å pixel −1 and a wavelength coverage of 4000 -7000Å. We adopted a slit width of 5 in order to minimize the light losses due to the aperture effect. We reduced the spectra with IRAF v2.16, and extracted them using a uniform aperture of ±6 .8 and background windows of 7 .6 -15 .1 on both sides. The wavelengths of the spectra are calibrated by taking CuAr lamp exposures. For each object, we took several consecutive exposures in order to estimate the systematic flux calibration uncertainties (see Section 2.2.2).

Spectrophotometry
We initially flux calibrated the spectra using one or more spectrophotometric standard stars observed each night (primarily Feige 34, G191B2B, or BD+28d4211). However, due to variable atmospheric extinction during the night, we took additional measures to obtain accurate photometry. We used established techniques to use [O III] lines for flux calibration (van Groningen & Wanders 1992; Fausnaugh 2017), which ensures good accuracy even in relatively poor observing conditions. The variable time scales of the [O III] λ5007 lines are much longer than one year for luminous AGNs as they originate from much more extended narrow-line regions (Peterson et al. 2013). Therefore, [O III] lines can be used as flux standards to calibrate the spectra.
We used a 5 -wide slit, which is wider than the FWHM of the seeing (2 ∼ 3 ) during all of the observations. The variable seeing at different times leads to the change of the line spread function (spectral resolution, see more details in Du et al. 2016a profiles by subtracting the local continuum underneath determined by interpola-tion between two nearby background windows. We provide the extraction and local continuum windows in Table 2. The optimal parameters of ζ were determined by the Levenberg-Marquardt algorithm. Then, all of the spectra were smoothed by convolution with their corresponding ζ to minimize the influence due to the variable seeing and the tracking inaccuracies. It should be noted that the spectral resolution after the convolution is lower than that of the original spectrum, but still quite sufficient to resolve broad Hβ profiles. After that, each exposure of the object was scaled to match its standard [O III] flux. The [O III] fluxes were measured using the windows listed in Table 2, and the standard [O III] fluxes (listed in Table 2) of the objects are determined by the spectra taken in the photometric conditions. We produced the final calibrated spectra by averaging the (appropriately noise-weighted) exposures in the same night for each object.

Mean and RMS spectra
To demonstrate the spectral characteristics and to evaluate the variation amplitude at different wavelengths, we plot the mean and root-mean-square (rms) spectra for each object in Figure 2. The mean and rms spectra are defined as and respectively. F i λ is the i-th spectrum, and N is the number of spectra for this object. The [O III] and the narrow Hβ emission lines in the rms spectra are extremely weak, which indicates that the calibration procedure in Section 2.2.2 works very well. The obvious broad Hβ lines in the rms spectra imply that their variations are significant.

LIGHT CURVES 3.1. Light curves from WIRO
The fluxes of the Hβ emission line can be measured by direct integration (e.g., Peterson et al. 1998;Kaspi et al. 2000;Bentz et al. 2009;Grier et al. 2012;Du et al. 2014;Fausnaugh et al. 2017) or by spectral fitting (e.g., Barth et al. 2013;Hu et al. 2015). The first method measures the flux by integrating the Hβ line after subtracting the local background determined by two continuum windows. The second method separates the emission lines from the continuum by multicomponent spectral fitting, and has been gradually adopted in recent years. Considering that (1) it is difficult to fit the very complex and asymmetric Hβ profiles of our targets perfectly by multiple Gaussian or Lorentzian functions, or their combinations; (2) the integration method is more robust and works well for isolated emission line like Hβ, we choose to use the integration method to measure the Hβ light curves in this work.
We chose the windows for Hβ flux measurements that cover the Hβ emission shown in their rms spectra ( Figure 2) and avoid the possible influences from their [O III] lines. The local continuum windows were selected as minimally variable regions in the rms spectra. We provide the line and the local continuum windows of Hβ in Table 2 and show them in Figure 2. We measured the Hβ light curves by integrating the fluxes in the Hβ windows after subtracting the local continuum determined by interpolating between the two nearby continuum windows. Similar to the narrow [O III] lines, the flux of the narrow Hβ can be regarded as a constant during our campaign. Thus, we did not remove the contributions of the narrow Hβ lines from the measured Hβ light curves. We obtained the 5100Å continuum light curves by averaging the fluxes from 5075 to 5125Å in the rest-frames (shown as grey regions in Figure 2).
The error bars of the fluxes in the light curves include two components: (1) the Possion noise and (2) the systematic uncertainties. The systematic uncertainties result primarily from the variable atmospheric extinction (especially in poor weather conditions) and telescope-tracking inaccuracies, and we estimated them using the scatter of the mean fluxes of the exposures in the same night over a wider range of wavelength (∼4700-5100Å). The two components summed in quadrature provide the error bars of the points in the light curves (see Figure 3). However, the above error bars likely do not account for all of the systematic uncertainties for most of the objects (the flux differences between adjacent nights are sometimes larger than the error bars and unlikely to result solely from real variability). We show additional systematic uncertainties estimated by the median filter method (see more details in Du et al. 2014), which smooths the light curve by a median filter of 5 points and then obtains the systematic uncertainty from the standard deviation of the residuals after subtracting the smoothed light curve, if necessary, as the gray error bars in the lower-right corners in Figure 3 3 . These are also taken into account in the following time-series analysis (Section 4.1) by added in quadrature to the error bar of every data point in Figure 3. We provide the 5100Å and Hβ light curves in Figure 3 and Table 3.

Photometric light curves from ASAS-SN
Photometric observations based on imaging are commonly carried out simultaneously with the spectroscopic observations in many RM campaigns (e.g., Bentz et al. 2009;Denney et al. 2010;Du et al. 2014Du et al. , 2015Du et al. , 2016bDu et al. , 2018Wang et al. 2014;Fausnaugh et al. 2017). The photometric light curves can be used to check the calibration precision of the spectroscopic observations, and furthermore can be adopted as sup-  plemental to the 5100Å light curves, especially if the sampling of the spectroscopic observations is relatively poor or their monitoring period is not long enough. We did not conduct our own photometric observations during 2016-2017 in our campaign. However, because the objects of this paper are bright enough, we can find their photometric light curves from archival data from the All-Sky Automated Survey for Supernovae (ASAS-SN) 4 . ASAS-SN is a long-term project useful to discover and study supernovae, transients, 4 http://www.astronomy.ohio-state.edu/ ∼ assassin/index.shtml and other variable sources by automatic and regular sky survey, and provides photometric light curves for the objects down to ∼17 magnitude. More information and technical details about the ASAS-SN light curves are provided by, e.g., Shappee et al. (2014) and Kochanek et al. (2017). Figure 3 shows the scaled ASAS-SN light curves of our targets (more details of the scaling are provided in Section 3.3). We removed several points with very poor S/N. The variations of the 5100Å and the ASAS-SN light curves are consistent (see Figure 3), thus it verifies our spectroscopic calibration.

Combined continuum light curves
Considering that the ASAS-SN observations can extend the temporal spans of our continuum light curves and improve their sampling cadences, we averaged the ASAS-SN and the 5100Å light curves to produce a combined continuum light curve for each object. Because of differing apertures (ASAS-SN uses an aperture with a radius of 15 .6), the ASAS-SN light curves require adjustment to match the mean fluxes and the variation amplitudes of the WIRO 5100Å light curves before combination. This was performed by assuming where F 5100 and F ASAS−SN are the 5100Å and ASAS-SN fluxes of the closely adjacent pairs of the observations (the separation is at most less than 2 days), a is a flux adjustment, and b is a scale factor. We determined the values of a and b by using the FITEXY algorithm (Press et al. 1992). Then, we scaled the ASAS-SN light curves by applying Equation √ " means we use this time lag of the object to calculate its black hole mass in Table 6. (5), and combined the 5100Å and ASAS-SN light curves by weighted averaging all of the observations in the same nights. The uncertainties of the weighted mean 5 are simply adopted as the error bars in the combined light curves. We tried to use the median filter method to estimate the systematic uncertainties, and found no extra systematic uncertainties are needed for the combined light curves. Figure 3 shows both the scaled ASAS-SN light curves and the final combined light curves.

ANALYSIS 4.1. Cross-correlation function
The time delays between the variations of the continuum and Hβ emission lines were determined by using the interpolated cross-correlation function (ICCF; Gaskell & Sparke 1986;Gaskell & Peterson 1987). To measure the time delay we used the centroid of the ICCF above 80% of the peak (a typical value used in many RM investigations, e.g., Bentz et al. 2009;Du et al. 2014;Fausnaugh et al. 2017).
We estimated the uncertainties of the time delays through the "flux randomization/random subset sampling (FR/RSS)" method (Peterson et al. 1998. The procedure takes into account both the measurement errors of the fluxes and the uncertainties due to the sampling/cadence. This method generates light curve realizations by perturbing the fluxes in accordance with their error bars and randomly sub-sampling the data points in the light curves. The cross-correlation centroid distributions (CCCDs) are obtained by performing the ICCF to the light curve realizations. The median and the 68% confidence intervals of the CCCDs are adopted as the final time lags and their uncertainties. The auto-correlation functions (ACFs), the CCFs, and the CCCDs of the 5100Å and Hβ light curves for each object are shown in Figure 3. Table 4 gives The time lag measurements and the correspond-5 The uncertainties of the weighted mean is defined as σmean = (1/ σ −2 i ) 1/2 , where σ i is the uncertainty of each point in the same night. ing maximum correlation coefficients of the CCFs. We also measured the Hβ time delays relative to the combined continuum light curves (see Section 3.3), and provide the results in Figure 3 and in Table 4. Although the photometric data can in principle extend the temporal spans of the 5100Å light curves and improve the sampling cadences, the scatter of the ASAS-SN light curves are generally larger than that of the 5100Å light curves in our campaign because of the limited collecting area of the ASAS-SN telescopes. We used the lag measurements of the combined continuum versus the Hβ light curves only in the case of Mrk 6, which greatly extended continuum coverage and resulted in a clearly improved result. The Hβ time lags we used in the following analysis are labeled by " √ " in Table  4.

Line-width measurements
We measured the line width (using both the FWHM and line dispersion σ Hβ parameters) of the broad Hβ for both the mean and rms spectra. The narrow-line contributions in the rms spectra are fairly weak for all of the targets, thus the corresponding FWHM and σ Hβ of the rms spectra can be obtained easily. We measured the Hβ profiles after subtracting continuum; continuum windows were selected beyond any contribution from the emission lines and can differ between the mean and rms profiles and from the ones listed in Table  2, which were optimized for light curves. However, before measuring the Hβ widths from the mean spectra, the narrow Hβ and [O III]λ4959, 5007 lines still do need to be removed.
We first assumed that all of the three narrow lines (Hβ and    present sample is very complex and asymmetric, we did not fit the entire profile but a narrower and local window (usually a window of 4000 ∼ 5000 km s −1 ) around 4861Å. We tried to add more Gaussians or changed the Gaussian(s) to a high-order polynomial as the broad Hβ contribution in the fitting, but the fitting results and the following width measurements did not change significantly because the shape of the narrow template has been constrained by the corresponding [O III]λ5007 profile of each object.  Stern & Laor 2013) in AGNs, which means that the narrow-line subtraction procedure adopted here appears to work well. Then we measured the FWHM and σ Hβ from the mean Hβ profiles after the narrow-line subtraction. We estimated the broad Hβ velocity width uncertainties for both the mean spectra and the rms spectra by using the bootstrap method. A subset spectrum was created by resampling N points with replacement from the N data points in the mean/rms spectrum. For the rms spectrum, we measured the FWHM and σ Hβ from the subset spectrum, and repeated the resampling and the measurement 500 times to generate the FWHM and σ Hβ distributions. The resulting median values and the standard deviations of the distributions were regarded as the measurements and the uncertainties. For the mean spectrum, we subtracted the narrow lines from the subset spectrum by using the procedure described above before the measuring the widths. The uncertainties of the narrowline flux ratios in Table 5 were also obtained at the same time.
To estimate the width of the line spread function (instrumental broadening) in our observations, we measured the widths of the [O III] lines in the mean spectra and compared them with the intrinsic narrow-line widths in Whittle (1992) or the higher resolution spectra of the Sloan Digital Sky Survey (SDSS). The FWHM of the line spread function obtained for different object ranges from ∼850 km s −1 to ∼1000 km s −1 . For simplicity, we adopted the mean value of 925 km/s (FWHM) as the line spread function for all of the objects in our campaign. After correcting the contribution of the line spread function, the line widths (FWHM and σ Hβ ) of the Hβ in the mean and rms spectra are listed in Table 6.

Black hole masses and accretion rates
Combining the time lag with the line width measured from the FWHM or line dispersion, the black hole mass can be obtained by the Equation 1. For AGNs with extreme BLR kinematics, e.g., super-fast or ultra-strong inflow or outflow, or extreme inclination angles, the virial factor f BLR may possibly differ significantly from typical values (e.g., Pancoast et al. 2014).
The mean f BLR of a sample can be calibrated by comparing the RM objects, which have bulge stellar velocity dispersion measurements (σ * ), with the M • -σ * relation of the inactive galaxies (e.g., Onken et al. 2004;Woo et al. 2010Woo et al. , 2015Park et al. 2012;Grier et al. 2013a;Ho & Kim 2014, see a brief review in Du et al. 2017). However, each individual object may have a very different virial factor (e.g., Horne 1994; Horne et al. 2004;Bentz et al. 2010;Grier et al. 2013b;Pancoast et al. 2014), especially for those AGNs with asymmetric Hβ profiles which may host BLRs with complex geometry or kinematics. We adopt f BLR = 1.12, 4.47 in Woo et al. (2015) corresponding to the FWHM and σ Hβ in the rms spectra, respectively, and also provide simple virial products (assuming f BLR = 1) for the FWHM measurements in the mean spectra for the present sample in Table 6, but acknowledge the potentially large uncertainty on f BLR . The uncertainties of the black hole masses listed in Table 6 only account for the error bars of the lag and width measurements.
We provide general estimates of their dimensionless accretion rates, defined asṀ =Ṁ • /L Edd c −2 , whereṀ • is the mass accretion rate and L Edd is the Eddington luminosity (Du et al. , 2016b. TheṀ can be estimated by using the formula (see more details in Du et al. 2015Du et al. , 2016b) where 44 = L 5100 /10 44 erg s −1 is the monochromatic luminosity at 5100Å, m 7 = M • /10 7 M , and i is the inclination angle of disk to the line of sight. We adopted cos i = 0.75 (see some discussions in Du et al. 2016b), which is an average estimate for type I AGNs (e.g., Pancoast et al. 2014). For the most precise results, it is necessary to subtract the host contribution from L 5100 before calculating the accretion rates, but this is beyond the present scope of this paper, so our estimates are upper limits. We foundṀ 1 ∼ 2 for 3C 120 and SBS 1518+593, and 0.2 for Ark 120 and Mrk 6. Therefore, all of the 4 objects are sub-Eddington accretors. More detailed determinations of luminosity will be done in a future paper allowing more precise estimates.

Velocity-resolved time lags
In order to investigate the geometry and kinematics of their BLRs, we measure the velocity-resolved time lags of the Hβ emission lines of our target AGNs. Several typical transfer function models and their corresponding velocityresolved time lags are given by, e.g., Bentz et al. (2009), and more complicated examples of the transfer function are provided by, e.g., Welsh & Horne (1991) and Horne et al. (2004). In general, longer lags in the high-velocity blue/red part of the emission line are regarded as the signatures of inflow/outflow, while a symmetric velocity-resolved structure around zero velocity, with smaller time lags for higher velocities, is diagnostic of Keplerian disk or at least virialized motion over a spatially extended BLR in general. We divided the Hβ lines into several velocity bins, each of which having the same integrated fluxes in their individual continuumsubtraced rms spectra based on interpolation between the windows shown in Figure 2. Then, we measured the light curve in each bin and performed an ICCF analysis (using the method in Section 4.1) with the 5100Å continuum light curve

NOTE-
The line spread function caused by the instrument and seeing has been corrected from the line-width measurements. MVP is the virial product measured from the mean spectrum (see more details in Section 4.3). The black hole masses (M•) are calculated from the rms spectra using the fBLR factors in Woo et al. (2015).
(the combined continuum light curve in the case of Mrk 6). Figure 4 shows the resulting time lags as a function of the velocity and the corresponding rms spectrum for each object. The velocity bins with equal rms flux have the same level of variation but may have different amounts of physical flux. As a further test, we divided the Hβ lines into the velocity bins, each of which having the same Hβ fluxes in the narrowline-subtracted mean spectra, and measured the velocityresolved lags in Appendix B. In general, the results are almost the same as the rms-based velocity-resolved lags, which means the velocity-resolved analysis here is robust. In the following Section 5, the discussion of their BLRs is based on the results in Figure 4.

DISCUSSION
Our data have produced integrated Hβ time lags as well as high-quality velocity-resolved time lags for our four featured AGNs. We discuss each object individually below and compare our results to past work as appropriate.

3C 120
As a nearby broad-line radio galaxy, its mean Hβ profile is asymmetric toward the red, however the rms Hβ profile is strongly blueshifted (see Figures 2 and 4). Its Hβ time lag respect to the varying continuum was first detected successfully by Peterson et al. (1998), albeit with large uncertainty, and the measured lag was 43.8 +27.7 −20.3 days in the rest frame. After ∼ 11 years, it was monitored again in 2008 -2009, and the observed rest-frame Hβ time lag was 27.9 +7.1 −5.9 days (Kollatschny et al. 2014). 3C 120 was observed the third time with higher temporal sampling in 2010 -2011, and a similar Hβ delay of 25.9 +2.3 −2.3 days in the rest frame was obtained from the light curves by using CCF analysis method . From the velocity-resolved time lag measurement and the transfer function reconstructed by the maximum entropy method (Grier et al. 2013b), its BLR was likely an inclined disk or a spherical shell, but there was also some evidence of inflow given that the strength of its line response was asymmetric (Grier et al. 2013b).
Our campaign was carried out ∼7 years later than the observation in Grier et al. (2012), and captured a very strong peak around Julian date ∼90 days (from the zero point of 2457700 in Figure 3) in the continuum and a clear response around ∼110 days in the Hβ light curve. The detected Hβ lag is 20.2 +5.0 −4.2 days in the rest frame, which is slightly shorter than the value in Grier et al. (2012), but the difference is not significant considering the uncertainties. In addition, the rms spectrum in our campaign is different from that in . Our rms spectrum is significantly blueshifted (see Figure 4), while the rms spectrum in Grier et al. (2012) has a strong red asymmetry. This implies that the varying part of the BLR in 3C 120 has significantly changed after ∼7 years.
Furthermore, the velocity-resolved lag measurement shows a complicated structure, which is different from the symmetric velocity-resolved results of a inclined disk or a spherical shell determined by Grier et al. (2013b). From 1500 km s −1 to −1500 km s −1 , the lag gradually decreases, which is the signature of outflow. However, the tendency changes around ∼ 1500 km s −1 , and the lags become longer at the blue end. This complicated structure suggests that its BLR is in a complex state. Of course, it should be noted that only the two bins (with a little larger uncertainties) at the highest blue velocities are in charge of the upturn at the blue side. More observations with better spectral resolution and higher S/N ratios are needed in the future to verify this complex state and to investigate the detailed BLR kinematics.

Ark 120
Given its asymmetric/double-peak Hβ profile and the longterm periodic-like variation in past decades, Ark 120 is a possible candidate for a supermassive black hole binary system (Li et al. 2017). Its Hβ time lags during 1990 September -1991March and 1995September -1996April, measured by Peterson et al. (1998, were 47.1 +8.3 −12.4 days and 37.1 +4.8 −5.4 days in the rest frame, respectively. Doroshenko et al. (2008) monitored this object from 1992 to 2005, and obtained an Hβ lag of 70 ± 7 days in the rest frame, by combining their data with other additional data from 1988 to 1996 (Peterson et al. 1998).
During our campaign, the mean spectrum of Ark 120 was not significantly double-peaked, but the rms spectrum shows clearly two peaks. One is located at roughly 4861Å in the rest frame, and the other is strongly redshifted by ∼ 2500  Table 4. km s −1 . The rms spectrum is different from the one in Peterson et al. (1998). The rms spectrum in Peterson et al. (1998) has three peaks (an additional blueshifted peak), and the redshifted peak is relatively stronger. One simple and reasonable guess is that both of the blueshifted and redshifted peaks became relatively weaker with respect to the peak with zero velocity by the epoch of our campaign, which implies a possible correlation between the origins of these two peaks. Our velocity-resolved lag measurement of Ark 120 is also complicated. In general, the time lag decreases from the blue (−3000 km s −1 ) to the red (4000 km s −1 ) velocity, which is the signature of an inflowing BLR. However, there is a local peak around 1000 ∼ 2000 km s −1 , which corresponds to the dip between the two peaks in the rms spectrum. To further investigate the geometry and kinematics of the BLR in Ark 120, reconstructing its velocity-delay map is needed. We will reconstruct the velocity-delay map of the objects in the present sample using the maximum entropy method (e.g., Horne 1994) in a separate paper in the near future.

Mrk 6
It has been known for more than 40 years that the Hβ profile of Mrk 6 is extremely asymmetric and has a blueshifted peak with the velocity of ∼ −3000 km s −1 (Khachikian & Weedman 1971). With its strange radio morphology (jet flips and jet precession as reported by, e.g., Kukula et al. 1996;Kharb et al. 2014), Mrk 6 is suggested to be a potential supermassive binary black hole system (Kharb et al. 2014).
Mrk 6 has been spectroscopically monitored by different groups in the past decades (Sergeev et al. 1999;Doroshenko & Sergeev 2003;Doroshenko et al. 2012;Grier et al. 2012). With the long temporal span of the campaign (1992 -2008), Doroshenko et al. (2012) found the time lag of its Hβ line to be 21.1 ± 1.9 days in the rest frame. Grier et al. (2012) observed this object again in 2010 with a higher sampling rate, and obtained a different Hβ lag of 10.1±1.1 days. In our observation, the new Hβ time lag is 18.5 +2.5 −2.4 days in the rest frame, which is more similar to the result of Doroshenko et al. (2012). Moreover, the line width (FWHM= 5274 km s −1 ) in our rms spectrum is similar to the value (5445 km s −1 ) in Doroshenko et al. (2012), and much smaller than the number (9744 km s −1 ) in Grier et al. (2012). These changes might be simply due to the BLR "breathing" effect (e.g., Peterson et al. 2002;Korista & Goad 2004;Cackett & Horne 2006).
Considering the very complex Hβ profile of this object, it is interesting to compare the rms spectrum and the velocityresolved time lags in our campaign with those reported by Doroshenko et al. (2012) and Grier et al. (2012). Doroshenko et al. (2012) found that the rms spectrum had two prominent peaks (one is blueshifted and the other has roughly zero velocity) during 1993 -1999. The two peaks in the rms spectrum almost disappeared in -2002, and rose again in 2005(Doroshenko et al. 2012. A third small peak with redshifted velocity (∼ 1500 km s −1 ) appeared in 2005 -2008, which makes the things even more complex. Grier et al. (2012) found that the blueshifted peak became strong again in 2010, and the redshifted peak was also still significant. It should be noted that the blueshifted peak always stayed at ∼ −3000 km s −1 , and did not show large velocity changes, however, its relative intensity has changed significantly with time (Doroshenko et al. 2012;Grier et al. 2012).
The velocity-resolved time lag measurements of Mrk 6 also shows very complex structures (Doroshenko et al. 2012;Grier et al. 2013b). Doroshenko et al. (2012) found that the velocity-resolved lags are generally shorter in the wings and longer in the line core, but the longest lag is blueshifted by ∼ −2000 km s −1 . They interpreted this complex velocityresolved lag measurement as a combination of virial motion plus inflowing gas. The velocity-resolved lags in Grier et al. (2013b) are similar to the results for 1993 -1995 in Doroshenko et al. (2012), but the lags increase gradually from 2000 km s −1 to 7500 km s −1 . Our velocity-resolved lags are more similar to the results of Grier et al. (2013b), that the longest lag is located at −2000 km s −1 and the lags increase from 1500 km s −1 to 7000 km s −1 . It is not easy to simply explain this complicated velocity-resolved lag structure using a combination of virial motion and inflowing gas suggested by Doroshenko et al. (2012). Thus, the velocitydelay map and detailed modeling of a more complex BLR geometry and kinematics are needed to further explore the nature of this object. It should be noted that the bin with the highest velocity on the red side has relatively larger uncertainty, which may be caused by the slight influence from the [O III]λ4959 line. We will conduct the observation in the future to check this issue.
Mrk 6 has the longest monitoring period, the highest sampling rate, and the best S/N ratios among the 4 objects in the present paper. Although the continuum of Mrk 6 changes very slowly (rise slowly in the first half and then shows a gently falling trend, see panel c in Figure 3), which makes its ACF and CCF exceptionally broad (see also Figure 3), the uncertainties of the its time lag measurement produced by the FR/RSS method (Peterson et al. 1998 are still acceptable. But it should be pointed out that the very slow variation in the continuum and the corresponding broad ACF/CCF may potentially limit the smallest observable lag, either probing different parts of the BLR or skewing the gas distribution to respond to larger radii (Goad & Korista 2015). The future observation will investigate this possibility.

SBS 1518+593
The time lag of SBS 1518+593 is mainly determined by the dip around 105 days in the continuum light curve and the response around 125 days in the Hβ light curve. The peak at ∼150 days in the continuum and its potential response at ∼180 in the emission-line light curve also provide constraint to the lag measurement, but only play a minor role because the peak in the Hβ light curve is already near to the end of the campaign.
The Hβ emission line in the mean spectrum of SBS 1518+593 is asymmetric toward the red wing, whereas the peak of the Hβ in the rms spectrum is mildly redshifted (see Figure 2). However, its velocity-resolved lag measurement is not significantly asymmetric. The longer lags in the line core and the shorter lags in the line wings imply that the geometry and kinematics of its BLR tend to be virialized motion or a Keplerian disk.

Ongoing project
Our primary objectives are: (1) revealing the complex BLR physics behind AGNs with asymmetric Hβ, (2) understanding the influence of differing BLR geometry or kinematics for the black hole mass measurement, and (3) looking for SMBH binary systems.
In particular, with the data quality and calibration precision improving in recent years, the velocity-resolved RM, including the velocity-resolved time lag measurements (e.g., Bentz et al. 2008Bentz et al. , 2009Bentz et al. , 2010Denney et al. 2009Denney et al. , 2010Grier et al. 2013b;Du et al. 2016a) and the analysis of the velocity-delay maps reconstructed through the maximum entropy method (e.g., Horne 1994;Horne et al. 2004;Bentz et al. 2010;Grier et al. 2013b), regularized linear inversion (Skielboe et al. 2015), or Bayesian-based dynamical modeling by Markov Chain Monte Carlo (MCMC) method (e.g., Pancoast et al. 2011Pancoast et al. , 2012Pancoast et al. , 2014Grier et al. 2017a), was successfully applied to more than 20 AGNs, and preliminarily revealed the geometry and kinematics of their BLRs. We plan to significantly add to this total and help make a breakthrough in our understanding.
6. SUMMARY In this paper, we describe the MAHA project and report some results from the first campaign. We successfully obtained the Hβ time lags of 4 objects observed from December 2016 to May 2017, and preliminarily investigated their BLR kinematics through measuring the velocity-resolved lags. The velocity-resolved results of 3C 120, Ark 120, and Mrk 6 showed very complex structures, which were different from the simple signatures of outflow, inflow, or virialized motion. The velocity-resolved lag measurements of SBS 1518+593 showed generally shorter lags in the line wings and longer ones at the line centers, which implied that their BLR is virialized motion or a Keplerian disk. The complexities of the velocity-resolved time lags in the AGNs with asymmetric Hβ line profiles clearly demonstrate the very complex geometry and kinematics of their BLRs, and provide good opportunities to understand the physics of the BLRs in AGNs in more details in the future.
We thank the anonymous referee for constructive comments. We acknowledge the support by National Key R&D Program of China (grants 2016YFA0400701 and 2016YFA0400702), by NSFC through grants NSFC-11503026, -11233003, -11573026, -11773029 Table 7, and plot recent WIRO spectra in Figure  5. The WIRO spectra displayed in Figure 5 have had their narrow lines subtracted as described in the main text.
Below are notes on individual objects in alphanumeric order by name. 2MASX J21140128+8204483. Also known as S5 2116+81, the Hβ profile displays a significant blue asymmetry. 3C 120. One of the objects featured in this paper and possessing an Hβ line with a red asymmetric tail. 3C 120 has been previously reverberation mapped several times, with recent results published by Grier et al. (2012).  3C 390.3. Others have previously reverberation mapped this object, which displays double-peaked emission lines (e.g., Shapovalova et al. 2010;Dietrich et al. 2012;Sergeev et al. 2017). 3C 390.3 once displayed quasi-periodic emission-line profile changes suggestive of a binary (Gaskell 1996), but later deviated from the predicted pattern (Eracleous et al. 1997).
III Zw 2. This object has an Hβ profile with an extreme red asymmetry. Ark 120. One of the objects featured in this paper, displaying a broad and complex double-peaked Hβ profile. Previous RM results exist in the literature (e.g., Peterson et al. 1998;Doroshenko et al. 2008).
Mrk 6. One of the objects featured in this paper and possessing a highly blueshifted peak and a red tail, probably also doublepeaked at the present epoch. There has been high-quality RM at previous epochs (Doroshenko et al. 2012;Grier et al. 2012).
Mrk 715. Also known as SDSS J100447.61+144645.6, this object has a double-peaked Hβ line profile and a long tail to the red.
Mrk 876. The current Hβ profile is suggestive of a double-peaked profile with red asymmetry. In the past (De Robertis 1985), the blueshifted peak was significantly stronger, reminiscent of Mrk 6.
Mrk 1148. This object has an Hβ profile with a mild red asymmetry. NGC 985. This object has an Hβ profile with a red asymmetry. PG 0947+396. The Hβ profile shows a red asymmetry. PG 1004+130. This object is somewhat luminous (> 10 45 ergs s −1 ) and a likely radio-loud broad absorption line quasar (Wills et al. 1999). The Hβ line has a red asymmetry.
PG 1048+342. This PG quasar has a red asymmetric Hβ profile. PG 1100+772. This luminous PG quasar has a Hβ profile with a blue bump but a red tail. PG 1151+117. The Hβ profile shows a red asymmetry. PG 1202+281. Also known as GQ COM, this object has an Hβ profile with a red asymmetry. PG 1302-102. This luminous radio-loud PG quasar has a red asymmetric Hβ profile. More notably, Graham et al. (2015) find a ∼5 year periodicity suggesting a binary nature, although Liu et al. (2018) suggest that the periodicity may have vanished. Time will tell. PG 1309+355. This PG quasar has a red asymmetric Hβ profile, and appears particular well fit by two Gaussians suggesting two components.
PG 1351+640. This object has an Hβ profile with a bump on the blue side but also a long red wing. The bump seems to have weakened compared to the spectrum shown by Boroson & Green (1992) observed over 25 years previously.
SBS 1518+593. One of the objects featured in this paper, featuring a mild red asymmetry at the present epoch. SDSS J015530.01-085704.0. The Hβ line has a significantly redshifted peak along with the customary associated blue asymmetry (Eracleous et al. 2012).
SDSS J023922.87-000119.5. The Hβ profile shows a red asymmetry. SDSS J093653.84+533126.8. In SDSS spectra, this object has an Hβ profile very similar to that of SDSS J094603.94+013823.6, showing a redshifted peak and blue asymmetry. Since then, a redshifted component has weakened dramatically leaving an emission line that is much more symmetric (Runnoe et al. 2017). All that is left now of that strong component is a weak, redshifted bump.
SDSS J094603.94+013823.6. The Hβ line has a significantly redshifted peak along with the customary associated blue asymmetry (Eracleous et al. 2012). The shifting profile is still consistent with expectations for a CB-SMBH (Runnoe et al. 2017).
SDSS J095539.81+453217.0. The Hβ line has a flat top with a blueshifted peak and a red asymmetry. SDSS J152139.66+033729.2. This object possesses an Hβ line with a red asymmetry. SDSS J171448.50+332738.2. This object possesses an unusual Hβ line with a redshifted top and red asymmetric wing. VIII Zw 233. The Hβ line has a significantly redshifted peak. WISE J134617.54+622045.3. The Hβ line has a significantly blueshifted peak along with the customary red asymmetry (Eracleous et al. 2012). The profile is reminiscent of Mrk 6. B. VELOCITY-RESOLVED TIME LAGS BASED ON MEAN SPECTRA The velocity-resolved lags in Section 4.4 are measured based on the velocity bins with equal flux in the rms spectrum, where the emission line in each bin have the same level of variation but not the physical flux. As a further test, we divided the Hβ lines into the velocity bins, each of which having the same Hβ fluxes in the narrow-line-subtracted mean spectra, and measured the velocity-resolved lags again. Figure 6 demonstrates the velocity-resolved lag measurements based on the mean spectra. Comparing with the rms-based results, the bins in high velocities become relatively narrower and the bins in low velocities become wider, because the Hβ profiles in the mean spectra are broader than those in the rms spectra for these objects (see Table  6). For Mrk 6, the red wing is located beneath the [O III]λ4959 (see Figure 6). We reduced the number of bins (compared with the rms-based result of Mrk 6) to make the bins wider that the highest velocity bin can exactly cover the [O III]λ4959 in order  Table 4. The black dashed lines in the lower panels are the narrow emission lines.
to avoid the potential influence from the strong [O III]. For the other 3 objects, the number of bins and the bluest and reddest velocity limits are the same as the rms-based results ( Figure 4). In general, the results are almost the same as the velocity-resolved lags in Section 4.4, which means the velocity-resolved analysis in this work is robust.