# Velocity-Resolved Ionization Mapping of Broad Line Region. I. Insights into
Diverse Geometry and Kinematics

Sha-Sha Li Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216,
Yunnan, People’s Republic of China  Key Laboratory for the Structure and
Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216,
Yunnan, People’s Republic of China  Center for Astronomical Mega-Science,
Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing
100012, People’s Republic of China  Key Laboratory of Radio Astronomy and
Technology, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District,
Beijing 100101, People’s Republic of China  Hai-Cheng Feng Yunnan
Observatories, Chinese Academy of Sciences, Kunming 650216, Yunnan, People’s
Republic of China  Key Laboratory for the Structure and Evolution of Celestial
Objects, Chinese Academy of Sciences, Kunming 650216, Yunnan, People’s
Republic of China  Center for Astronomical Mega-Science, Chinese Academy of
Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, People’s Republic
of China  Key Laboratory of Radio Astronomy and Technology, Chinese Academy of
Sciences, 20A Datun Road, Chaoyang District, Beijing 100101, People’s Republic
of China  H. T. Liu Yunnan Observatories, Chinese Academy of Sciences, Kunming
650216, Yunnan, People’s Republic of China  Key Laboratory for the Structure
and Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming
650216, Yunnan, People’s Republic of China  Center for Astronomical Mega-
Science, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District,
Beijing 100012, People’s Republic of China  J. M. Bai  Yunnan Observatories,
Chinese Academy of Sciences, Kunming 650216, Yunnan, People’s Republic of
China  Key Laboratory for the Structure and Evolution of Celestial Objects,
Chinese Academy of Sciences, Kunming 650216, Yunnan, People’s Republic of
China  Center for Astronomical Mega-Science, Chinese Academy of Sciences, 20A
Datun Road, Chaoyang District, Beijing 100012, People’s Republic of China  Key
Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, 20A
Datun Road, Chaoyang District, Beijing 100101, People’s Republic of China
Xiang Ji Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80
Nandan Road, Shanghai 200030, People’s Republic of China  Yu-Xuan Pang
Department of Astronomy, School of Physic, Peking University, Beijing 100871,
People’s Republic of China  Kavli Institute for Astronomy and Astrophysics,
Peking University, Beijing 100871, People’s Republic of China  Cheng Cheng
Chinese Academy of Sciences South America Center for Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101,
People’s Republic of China  CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101,
People’s Republic of China  Kai-Xing Lu Yunnan Observatories, Chinese Academy
of Sciences, Kunming 650216, Yunnan, People’s Republic of China  Key
Laboratory for the Structure and Evolution of Celestial Objects, Chinese
Academy of Sciences, Kunming 650216, Yunnan, People’s Republic of China
Center for Astronomical Mega-Science, Chinese Academy of Sciences, 20A Datun
Road, Chaoyang District, Beijing 100012, People’s Republic of China  Jian-Guo
Wang Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216,
Yunnan, People’s Republic of China  Key Laboratory for the Structure and
Evolution of Celestial Objects, Chinese Academy of Sciences, Kunming 650216,
Yunnan, People’s Republic of China  Center for Astronomical Mega-Science,
Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing
100012, People’s Republic of China  Rui Li School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou 450001, People’s Republic of
China Hai-Cheng Feng & H. T. Liu<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Broad emission lines of active galactic nuclei (AGNs) originate from the
broad-line region (BLR), consisting of dense gas clouds in orbit around an
accreting supermassive black hole. Understanding the geometry and kinematics
of the region is crucial for gaining insights into the physics and evolution
of AGNs. Conventional velocity-resolved reverberation mapping may face
challenges in disentangling the degeneracy between intricate motion and
geometry of this region. To address this challenge, new key constraints are
required. Here, we report the discovery of an asymmetric BLR using a novel
technique: velocity-resolved ionization mapping, which can map the distance of
emitting gas clouds by measuring Hydrogen line ratios at different velocities.
By analyzing spectroscopic monitoring data, we find that the Balmer decrement
is anticorrelated with the continuum and correlated with the lags across broad
emission line velocities. Some line ratio profiles deviate from the
expectations for a symmetrically virialized BLR, suggesting that the red-
shifted and blue-shifted gas clouds may not be equidistant from the
supermassive black hole (SMBH). This asymmetric geometry might represent a
formation imprint, provide new perspectives on the evolution of AGNs, and
influence SMBH mass measurements.

Active galactic nuclei (16), Time domain astronomy (2109), Reverberation
mapping (2019), Supermassive black holes (1663), Photoionization(2060)

††facilities: Lijiang: 2.4 m††software: DASpec (https://github.com/PuDu-
Astro/DASpec), PyCALI (Li et al., 2014), JAVELIN (Zu et al., 2011).

##  1 Introduction

Active Galactic Nuclei (AGNs) are fascinating astrophysical objects that
harbor supermassive black holes (SMBHs) at their centers. Accretion of matter
onto these black holes generates intense radiation across the electromagnetic
spectrum (Salpeter, 1964). Among the multi-component structure of AGNs, the
broad-line region (BLR) is a crucial component responsible for producing broad
emission lines due to the high-velocity motion of gas clouds near the central
black hole. Comprehending the geometry and kinematics of BLRs is essential not
only for probing the physical mechanisms associated with accretion, such as
inflows that fuel AGNs (Grier et al., 2013; Zhou et al., 2019) or outflows
related to disk winds (Emmering et al., 1992; Nicastro, 2000; Elitzur et al.,
2014), but also for accurately measuring the masses of SMBHs, $M_{\rm{BH}}$
(Pancoast et al., 2014; Mejía-Restrepo et al., 2018). Accurate values of
$M_{\rm{BH}}$ are vital for studying formation and evolution of galaxy, as
well co-evolution of SMBHs and host galaxies (Kormendy & Ho, 2013).

In the current picture of the AGN unification scheme, the BLR is modeled as a
clumpy, extended structure located between the accretion disk and dust torus
(Antonucci, 1993; Netzer, 2015). Regardless of whether the structure is clumpy
or smooth (Laor et al., 2006), the extended BLR geometry has been supported by
a variety of observational evidence. Reverberation mapping (RM) observations,
which measure the radius of the BLR based on the delayed response of emission
lines to continuum variations (Blandford & McKee, 1982), reveal the “breathing
effect” (Gilbert & Peterson, 2003; Goad et al., 2004; Cackett & Horne, 2006;
Wang et al., 2020), where the luminosity-weighted radius changes in response
to continuum variations. This is further supported by radial stratification of
the BLR as measured by multiline RM (Clavel et al., 1991; Bentz et al., 2010;
Feng et al., 2021b; Lu et al., 2022). More direct evidence comes from GRAVITY
phase interferometry which shows the line-emitting gas occupying the area from
the center out to the hot dust radius in three objects (Gravity Collaboration
et al., 2018; GRAVITY Collaboration et al., 2020, 2021). Although the overall
scale of the BLR has been unambiguously confirmed, the details of geometry and
kinematics remain an ongoing debate.

A number of velocity-resolved RM studies prefer symmetrical configurations of
the BLR, either disk-like or spherical. Shorter delays in the line wings are
typically expected from a virialized flow, while longer responses in the line
wings are often attributed to the presence of additional outflows or inflows
(Welsh & Horne, 1991; Bentz et al., 2009; Grier et al., 2013; Waters et al.,
2016; Lu et al., 2019; Bao et al., 2022; Villafaña et al., 2022).
Nevertheless, asymmetrical geometric models are widely adopted when fitting
the complex broad line profile. Chen & Halpern (1989) and Chen et al. (1989)
demonstrated that the double-peaked broad emission lines could be well-fitted
by a relativistic accretion disk. Subsequently, some non-axisymmetric
accretion disk models, such as an elliptical disk (Eracleous et al., 1995) or
perturbations in the disk (Newman et al., 1997; Gilbert et al., 1999; Storchi-
Bergmann et al., 2003), were proposed to explain the stronger red peak.
Moreover, accretion disk winds (Flohic et al., 2012), bowl-shaped BLRs (Goad
et al., 2012), and BLRs illuminated by anisotropic continuum source (Wanders
et al., 1995; Goad & Wanders, 1996) can also yield double-peaked profiles,
while theoretical calculations indicate that a spiral-arm geometry can
generate the velocity-delay maps obtained in RM observations (Du & Wang,
2023).

The broad hydrogen Balmer decrements (BDs) are widely employed to probe the
internal reddening of AGNs (Reynolds et al., 1997; Dong et al., 2008), where
various values can be associated with different degrees of dust extinction.
This internal reddening can influence the BDs by absorbing and scattering
light, thereby altering the observed BDs. On the other hand, photoionization
model simulations predict that the intrinsic BD should cover a range of values
within a parameter grid (Korista & Goad, 2004; Guo et al., 2020; Wu et al.,
2023). These values are dependent on physical conditions like gas density,
ionization parameter, and the spectral energy distribution of the ionizing
continuum. Moreover, considering that higher-order lines are theoretically
expected to originate at smaller radii on average (Korista & Goad, 2004), the
BD may undergo changes due to reverberation effects within the spatially
extended BLR.

Unfortunately, the actual dust and gas environments may vary substantially
between individual AGNs, making it extremely challenging to differentiate the
influences of these two physical mechanisms on the BD using statistical
samples of AGNs. However, the short-term (weeks to months) variability of BD
in single objects could potentially test the influence of ionizing photon
intensity on intrinsic BD, considering that changes in the obscuring materials
and gaseous state have much longer timescales. As the incident continuum flux
declines in a power-law manner with distance for an isotropic BLR, it becomes
possible to deduce a distance-dependent gas distribution through the velocity-
resolved BD. This technique essentially breaks the degeneracy between geometry
and kinematics according to the ionization distribution of the gas, so we
named it “velocity-resolved ionization mapping (IM)”. Combining velocity-
resolved IM and RM results can shed more light on the geometry and kinematics
of the BLR.

Data from multiline RM campaigns, typically collected over a several-month
period of spectroscopic monitoring, are ideal for the aforementioned purposes.
Here, we present a joint analysis of lags and BDs as a function of line-of-
sight velocity to deduce the geometry of the BLR. The paper is structured as
follows. The sample is introduced in Section 2. In Section 3, we present
analysis processes, including spectral fitting, time series analysis,
measurement of emission line flux ratios, and transfer function. The results
are given in Section 4. Discussions of the results are provided in Section 5.
Section 6 is a summary. In the Appendix, we examine the impact of
reverberation effects on velocity-resolved BD and the variability of dust
extinction.

Figure 1: Light curves, CCF results, and correlation between continuum and
emission line flux ratios for three targets. Arp 151, Mrk 50, and NGC 2617
belong to LAMP2008, LAMP2011, and Lijiang campaigns, respectively. The top
panels show light curves of the continuum, H$\alpha$, H$\beta$, and flux
ratios of H$\alpha$ and H$\beta$ ($F_{\rm H\alpha}/F_{\rm H\beta}$), and CCF
results. Black points denote the observed data, while the orange-shaded
regions illustrate the results reconstructed using JAVELIN. The $F_{\rm
H\alpha}/F_{\rm H\beta}$ values, derived from both observed data and JAVELIN
reconstructed results, are depicted by black points and orange bands,
respectively. The ACFs of the continuum are also present. The black lines
denote CCFs between the continuum and H$\alpha$, H$\beta$, and $F_{\rm
H\alpha}/F_{\rm H\beta}$, with their corresponding CCCDs portrayed by blue
lines. Vertical grey dashed lines represent instances where time lags are 0
days. The bottom panels show the relationship between the continuum and
$F_{\rm H\alpha}/F_{\rm H\beta}$. The $F_{\rm H\alpha}/F_{\rm H\beta}$ values,
obtained from the direct observations (depicted as black points), JAVELIN
reconstruction (represented as orange points), and time-lag-corrected JAVELIN
reconstruction (shown as green points), exhibit an anti-correlation with the
continuum.

##  2 Sample

We selected AGN samples with RM observations to investigate the flux ratios of
H$\alpha$ and H$\beta$ emission lines, their RM results, and their joint
analysis. However, simultaneous measurement of time delays for both H$\alpha$
and H$\beta$ emission lines has been realised by only a handful of AGNs
observations. Fewer than 10 AGN have simultaneous velocity-resolved RM
measurements for both H$\alpha$ and H$\beta$. Here we focus on the Lick AGN
Monitoring Projects in 2008 and 2011 (LAMP2008 (Bentz et al., 2009) and
LAMP2011 (Barth et al., 2015)), along with our recently observed changing-look
AGN (CL-AGN) samples from the Lijiang RM spectroscopic monitoring project
campaign for our research (Feng et al., 2021a, b; Li et al., 2022).

The LAMP2008 campaign, which was designed to measure $M_{\rm{BH}}$ for 13
nearby Seyfert galaxies, spanned over 64 nights from March to June 2008 using
the 3-meter telescope at Lick Observatory (Bentz et al., 2009). As this
project targeted low-redshift AGNs and the spectral data are publicly
available, we are able to perform a reanalysis of the broad H$\alpha$ and
H$\beta$ emission lines. The campaign also included photometric observations
in the $B$ and $V$ bands, which showed no significant time lags between them
(Walsh et al., 2009). To improve the photometry sampling, we utilized
PyCALI111https://github.com/LiyrAstroph/PyCALI, a Python package for
intercalibrating light curves in the $B$ and $V$ bands, which applies the
damped random walk (DRW) process to model observed variations of AGNs (Kelly
et al., 2009; Li et al., 2014). We then combined the calibrated results within
0.1-day intervals to represent the AGN continuum, as shown in Figure 1.

The LAMP2011 campaign, a follow-up project utilizing the same telescope, was
conducted from March to June 2011 at Lick Observatory. This campaign was
focused on the study of 15 low-redshift Seyfert 1 galaxies (Barth et al.,
2015). The sample with an apparent magnitude of $V\leq 17$ were chosen to
ensure high-quality spectra and to determine velocity-resolved RM results. For
the analysis of LAMP2011, we directly used the light curves provided in Table
4 of Barth et al. (2015), as the public spectra were unavailable.

Additionally, we obtained data through our CL-AGN program from 2018 to 2021.
This program was intended to investigate the BLR nature of CL-AGNs, which
included NGC 3516, NGC 2617, and NGC 4151 (Feng et al., 2021a, b; Li et al.,
2022). These observations were made using the 2.4-meter telescope at Lijiang
Observatory. These low-redshift samples allowed us to simultaneously detect
H$\alpha$ and H$\beta$ emission lines. In summary, the selected samples with
low redshifts and high spectral resolution enable us to reveal the RM results
associated with the broad emission lines.

Figure 2: Spectral fitting results. The upper-left panel shows the initial
fitting result for an individual spectrum. The upper-right panel shows the
second round fitting result for the mean spectrum of Arp 151. In both cases,
the black lines are the original spectra; red lines are the best-fitting
results; grey lines represent the host galaxy template; blue lines depict the
AGN power law; green lines represent the Fe ii template; magenta lines
represent the broad Balmer lines; orange lines represent the narrow lines;
cyan lines represent broad He i and He ii; the wavelength ranges appearing
between two dashed vertical lines represent regions that were not considered
during the fitting process. The lower panels display the fit residuals.

##  3 Analysis

###  3.1 Spectral Fitting

At optical wavelengths broad emission lines in AGN spectra are blended with
other components like the host galaxy, AGN power-law continuum, Paschen
continuum, narrow emission lines, and optical Fe ii emission. Therefore, to
extract information on these broad emission lines, we have applied multi-
component spectral fitting techniques, utilizing
DASpec222https://github.com/PuDu-Astro/DASpec, a tool based on the Levenberg-
Marquardt algorithm. For the CL-AGN sample, their spectra were flux calibrated
using comparison stars, with details of data reductions provided in Feng et
al. (2021a, b); Li et al. (2022). We adopted previous fitting models for these
objects (Feng et al., 2021a, b; Li et al., 2022), and the details of the
models have been omitted from this text.

For the LAMP2008 sample, both reduced and scaled spectra are publicly
available. The scaled spectra were calibrated according to the method
described in van Groningen & Wanders (1992), which relies on the flux of [O
iii] $\lambda$5007 measured from the reduced spectra. However, Hu et al.
(2016) provided a more accurate flux calibration through spectral fitting to
measure the same [O iii] $\lambda$5007 fluxes from reduced spectra. Therefore,
we recalibrated the LAMP2008 spectra following their method, conducting two
fitting rounds to obtain the final broad-line information. Before fitting, we
corrected the spectra for Galactic extinction and redshift using parameters
from the NASA/IPAC Extragalactic Database and the $R_{\rm v}=3.1$ extinction
law (Cardelli et al., 1989; O’Donnell, 1994).

In the first fitting round, we selected narrow wavelength windows spanning
4430-4600 Å and 4750-5550 Å, encompassing H$\beta$ but excluding He ii
$\lambda$4686\. Our model consisted of: (1) a power law representing the AGN
continuum (with a minor contribution from the Paschen continuum), (2) an
optical Fe ii template from Boroson & Green (1992), (3) a simple Stellar
Population model from Bruzual & Charlot (2003) to fit the host galaxy,
assuming solar metallicity and an age of 11 billion years, (4) a single
Gaussian for narrow lines like [O iii] $\lambda\lambda$4959,5007, and the
narrow H$\beta$ line, (5) a fourth-order Gauss-Hermite function for broad
H$\beta$. Additionally, we fixed the widths and shifts of narrow lines to [O
iii] $\lambda$5007 and set the [O iii] $\lambda\lambda$5007,4959 flux ratio at
3. We first used the above model to fit spectra observed on April 10, 2008
(with steady photometric night as mentioned in Bentz et al. 2009), deriving [O
iii] $\lambda$5007 fluxes (an example is shown in left panels of Figure 2).
The ratio of this measurement to the value listed in Table 3 of Bentz et al.
(2009) (after correcting for extinction) was used as the calibration factor.
Finally, these factors were applied to the reduced spectra for calibration.

In the second fitting round, we began by fitting the calibrated spectra.
Before fitting, we selected spectra with a signal-to-noise ratio (S/N) $\geq$
25 pixel-1 around 5100 Å, which were used for analysis and to create the mean
spectrum. The fitting windows covered the ranges of 4305-5750 Å and 5800-6800
Å, excluding the region around 5800 Å contaminated by Na i D residual noise
(as mentioned in Bentz et al. 2010). Within these windows, we could fit nearly
all broad emission lines from H$\gamma$ to H$\alpha$ and determine their
characteristics. The templates of the host galaxy and Fe ii were similar to
those used in the initial fitting. We also used fourth-order Gauss-Hermite
functions to model broad Balmer lines, including H$\alpha$, H$\beta$, and
H$\gamma$ (for NGC 5548, three Gaussians were separately used to fit them),
but the profile of broad H$\gamma$ was bound with broad H$\beta$. Broad He ii
and He i, as well as narrow lines, were individually modeled by a single
Gaussian. The widths and shifts of narrow lines with wavelengths below 6000 Å,
e.g., narrow H$\beta$, are tied to [O iii] $\lambda$5007, while those of
narrow emission lines with wavelengths above 6000 Å are tied to [O i]
$\lambda$6300 (or [S ii] $\lambda$ 6718 if [O i] $\lambda$6300 is weak). Also,
the flux ratios of [O iii] $\lambda\lambda$5007,4959 and [N ii]
$\lambda\lambda$6583,6548 were separately fixed to 3.0 and 2.96. We first
utilized these models to fit the mean spectrum (shown in right panels of
Figure 2) and obtained the widths and shifts of broad He ii and He i lines, as
well as the flux ratios of narrow lines relative to [O iii] $\lambda$5007 and
[O i] $\lambda$6300 (or [S ii] $\lambda$ 6718), respectively. Then, the
profiles of broad He ii and He i, and these narrow-line flux ratios measured
from the mean spectrum were fixed to fit individual spectra. Lastly, we
measured the fluxes and widths of each emission line in all spectra.

Based on [O iii] $\lambda$5007 widths, we convolved each spectrum to match the
maximum [O iii] $\lambda$5007 width to alleviate the impact of different
spectral spreads. We then applied the aforementioned fitting method to the
convolved spectra. After the fitting process, the broad emission line profiles
extracted from these convolved spectra were used in the subsequent analysis.

When analyzing the light curves, we adopted fitting errors for the LAMP2008
and Lijiang data. For LAMP2011, we directly used the errors provided in Barth
et al. (2015). We also added a systematic error to all light curves following
the method in Li et al. (2021). Using the fluxes and errors, we computed the
relative errors for each light curve, with mean values listed in Table 1. To
assess the characteristic variability of the light curves, we computed the
fractional variability amplitude ($F_{\rm var}$), as governed by the
methodology outlined in Rodríguez-Pascual et al. (1997):

$F_{\rm var}=\frac{(S^{2}-\triangle^{2})^{1/2}}{\langle F\rangle},$ (1)

where $S^{2}$ is the flux variance, $\triangle^{2}$ is mean squared error, and
$\langle F\rangle$ is the mean flux over the light curve. The $F_{\rm var}$
uncertainty is given by Edelson et al. (2002):

$\sigma_{\rm var}=\frac{1}{F_{\rm
var}}\bigg{(}\frac{1}{2N}\bigg{)}^{1/2}\frac{S^{2}}{{\langle F\rangle}^{2}}.$
(2)

where $\it N$ is the total number of flux measurements within the light curve.
The computed results for each light curve are compiled in Table 1.

Table 1: Measurements of variability amplitude and relative error

|  $F_{\rm var}$ (%) | Relative error (%) |
---|---|---|---
Object | ${\rm con}$ | ${\rm H\alpha}$ | ${\rm H\beta}$ | ${\rm con}$ | ${\rm H\alpha}$ | ${\rm H\beta}$ | $\rm\Delta BD$
LAMP2008 campaign
Mrk 142 | $5.01\pm 0.84$ | $4.47\pm 0.62$ | $2.46\pm 0.27$ | 1.34 | 4.05 | 2.65 | 0.03
SBS 1116+583A | $6.30\pm 0.83$ | $8.78\pm 1.31$ | $9.96\pm 0.98$ | 2.92 | 2.87 | 5.33 | 0.27
Arp 151 | $12.76\pm 1.51$ | $19.60\pm 2.17$ | $15.98\pm 1.38$ | 1.44 | 3.49 | 1.29 | 0.78
Mrk 1310 | $7.78\pm 1.17$ | $15.14\pm 1.88$ | $10.97\pm 1.04$ | 2.85 | 5.20 | 6.78 | 0.77
Mrk 202 | $12.18\pm 1.35$ | $16.97\pm 1.85$ | $3.77\pm 0.41$ | 1.60 | 3.03 | 3.58 | 0.33
NGC 4253a | $2.93\pm 0.63$ | $3.58\pm 0.72$ | $2.85\pm 0.31$ | 1.16 | 3.10 | 3.52 | 0.26
NGC 4748 | $8.82\pm 1.04$ | $6.87\pm 0.80$ | $4.98\pm 0.56$ | 2.31 | 2.68 | 1.93 | 0.11
IC 4218a | $7.40\pm 1.15$ | $21.16\pm 3.27$ | $7.23\pm 0.84$ | 3.77 | 4.04 | 12.09 | 1.96
MCG -06-30-15a | $5.18\pm 0.76$ | $4.75\pm 2.01$ | $3.37\pm 0.40$ | 1.72 | 2.94 | 7.99 | 0.09
NGC 5548 | $6.15\pm 0.77$ | $8.53\pm 1.12$ | $8.95\pm 0.80$ | 1.12 | 2.82 | 4.50 | 0.20
Mrk 290a | $3.44\pm 0.55$ | $3.97\pm 0.52$ | $3.61\pm 0.39$ | 1.33 | 2.60 | 2.13 | 0.12
IC 1198a | $3.17\pm 0.74$ | $6.81\pm 1.72$ | $3.27\pm 0.44$ | 2.46 | 3.46 | 8.05 | 0.21
NGC 6814 | $5.21\pm 0.74$ | $9.14\pm 1.18$ | $16.97\pm 1.81$ | 3.18 | 2.97 | 4.07 | 0.39
LAMP2011 campaign
Mrk 40 (Arp 151) | $13.99\pm 1.61$ | $19.60\pm 2.22$ | $29.14\pm 3.43$ | 6.25 | 1.97 | 0.82 | 0.91
Mrk 50 | $10.39\pm 1.24$ | $19.76\pm 1.90$ | $33.91\pm 3.34$ | 7.12 | 3.15 | 1.88 | 1.73
Mrk 279 | $3.77\pm 0.66$ | $7.22\pm 0.94$ | $10.88\pm 1.77$ | 6.46 | 2.50 | 1.97 | 0.56
Mrk 704 | $6.39\pm 0.85$ | $4.58\pm 0.58$ | $11.64\pm 1.73$ | 6.47 | 2.51 | 1.47 | 0.30
Mrk 1511 | $7.33\pm 1.34$ | $11.67\pm 1.33$ | $18.40\pm 2.28$ | 5.91 | 5.82 | 1.02 | 0.41
NGC 4593 | $15.01\pm 1.80$ | $22.68\pm 2.48$ | $36.89\pm 4.24$ | 9.80 | 5.16 | 0.65 | 1.62
Zw 229-015 | $12.23\pm 1.94$ | $24.07\pm 3.43$ | $24.20\pm 3.39$ | 6.82 | 4.33 | 7.57 | 1.05
Lijiang CL-AGN campaign
NGC 2617 | $7.29\pm 0.67$ | $15.27\pm 1.40$ | $4.45\pm 0.39$ | 1.25 | 1.43 | 2.81 | 1.79
NGC 3516 | $26.33\pm 2.31$ | $40.80\pm 3.80$ | $2.59\pm 0.27$ | 1.00 | 4.37 | 15.20 | 3.25
NGC 4151 | $3.49\pm 0.42$ | $5.30\pm 0.62$ | $7.53\pm 0.52$ | 1.59 | 1.09 | 1.38 | 0.17

Note. — The objects analyzed in this study were selected from the LAMP2008,
LAMP2011, and Lijiang campaigns. The LAMP2011 datasets incorporate
contributions from narrow lines.

a Those five objects, whose variability characteristics prevent lag
measurements from being conducted, are also illustrated in Figure 6 of Bentz
et al. (2009).

Table 2: Cross-correlation measurements

Object | $r_{\rm max,H\alpha}$ | $\tau_{\rm H\alpha}\rm(days)$ | $r_{\rm max,H\beta}$ | $\tau_{\rm H\beta}\rm(days)$
---|---|---|---|---
LAMP2008 campaign
Mrk 142 | 0.54 | $2.87_{-1.72}^{+7.49}$ | 0.58 | $9.44_{-5.19}^{+2.35}$
SBS 1116+583A | 0.58 | $3.33_{-0.74}^{+1.58}$ | 0.70 | $1.85_{-0.57}^{+1.50}$
Arp 151 | 0.93 | $7.30_{-1.12}^{+0.97}$ | 0.98 | $3.10_{-0.54}^{+0.76}$
Mrk 1310 | 0.84 | $3.76_{-0.57}^{+0.95}$ | 0.78 | $2.98_{-1.02}^{+0.71}$
Mrk 202 | 0.75 | $16.73_{-6.74}^{+2.98}$ | 0.76 | $10.69_{-2.54}^{+3.29}$
NGC 4253 | 0.28 | $6.19_{-4.86}^{+1.21}$ | 0.66 | $5.59_{-1.71}^{+1.53}$
NGC 4748 | 0.78 | $10.16_{-2.59}^{+0.99}$ | 0.79 | $6.08_{-1.73}^{+1.52}$
IC 4218 | 0.37 | $-1.10_{-4.97}^{+4.14}$ | 0.48 | $2.71_{-4.72}^{+2.58}$
MCG -06-30-15 | 0.48 | $-17.75_{-0.81}^{+8.36}$ | 0.39 | $-0.68_{-6.32}^{+6.99}$
NGC 5548 | 0.60 | $10.52_{-1.90}^{+1.31}$ | 0.68 | $5.63_{-1.62}^{+1.36}$
Mrk 290 | 0.42 | $12.83_{-18.17}^{-0.06}$ | 0.72 | $3.32_{-2.19}^{+3.28}$
IC 1198 | 0.39 | $-19.11_{+0.14}^{+20.82}$ | 0.26 | $0.88_{-9.34}^{+6.43}$
NGC 6814 | 0.67 | $9.28_{-1.64}^{+1.08}$ | 0.77 | $3.96_{-0.77}^{+0.90}$
LAMP2011 campaign
Mrk 40 (Arp 151) | 0.83 | $8.45_{-1.68}^{+1.91}$ | 0.94 | $5.26_{-0.92}^{+0.77}$
Mrk 50 | 0.93 | $14.66_{-1.40}^{+2.08}$ | 0.94 | $8.09_{-1.34}^{+1.05}$
Mrk 279 | 0.73 | $17.53_{-5.29}^{+2.41}$ | 0.83 | $16.85_{-2.66}^{+1.50}$
Mrk 704 | 0.73 | $3.15_{-1.27}^{+3.18}$ | 0.89 | $-0.20_{-1.95}^{+1.16}$
Mrk 1511 | 0.77 | $7.30_{-1.93}^{+2.90}$ | 0.84 | $5.99_{-0.75}^{+1.49}$
NGC 4593 | 0.67 | $6.54_{-2.12}^{+3.72}$ | 0.88 | $3.19_{-2.10}^{+1.69}$
Zw 229-015 | 0.86 | $16.24_{-4.60}^{+5.41}$ | 0.89 | $3.81_{-0.96}^{+3.37}$
Lijiang CL-AGN campaign
NGC 2617 | 0.68 | $7.88_{-0.93}^{+0.89}$ | 0.76 | $5.82_{-0.88}^{+1.00}$
NGC 3516 | 0.81 | $19.89_{-7.85}^{+2.01}$ | 0.85 | $7.36_{-1.59}^{+3.73}$
NGC 4151 | 0.81 | $5.00_{-3.77}^{+0.89}$ | 0.88 | $4.95_{-1.20}^{+2.10}$

Note. — The objects analyzed in this study were selected from the LAMP2008,
LAMP2011, and Lijiang campaigns. Using the CCF method, we calculated the
maximum cross-correlation coefficients ($r_{\rm max}$) and the corresponding
time lags for H$\alpha$ and H$\beta$ relative to the continuum. These values
are denoted as $r_{\rm max,H\alpha}$, $r_{\rm max,H\beta}$, $\tau_{\rm
H\alpha}$, and $\tau_{\rm H\beta}$, respectively.

Figure 3: Correlation between velocity-resolved emission line flux ratios and
time lags of Arp 151 (from LAMP2008) and NGC 2617. In the top panels, the
black step lines represent the rms spectra of H$\alpha$. The rms spectra are
divided into multiple velocity bins, ensuring a uniform integrated rms flux
across each bin. Vertical dashed lines delineate the boundaries of these
velocity bins. The $F_{\rm H\alpha}/F_{\rm H\beta}$ values and time lags
($\tau_{\rm cent}$) of H$\alpha$ for each velocity bin are shown as black
points, with the velocity associated with each point aligning with the
centroid of the rms spectra in its respective bin. Each velocity bin was not
corrected for the delay of emission line relative to the continuum before
choosing the flux ratio. Horizontal dotted lines and grey bands depict the
time lags of H$\alpha$ and their associated uncertainties, respectively. In
the bottom panels, the relationship between $F_{\rm H\alpha}/F_{\rm H\beta}$
values and time lags of H$\alpha$, as obtained from various velocity bins.

Figure 4: Same as Figure 3, but for SBS 1116+583A, NGC 3516, and NGC 4151.

###  3.2 Time Series Analysis

The light curves of the sample allow us to estimate the time delays between
broad emission lines and the AGN continuum using the interpolated CCF method
(Gaskell & Sparke, 1986; Gaskell & Peterson, 1987). These time delays (denoted
as $\tau_{\rm cent}$) are determined by measuring the centroid of the segment
above 80% of the peak value in the CCF, with the peak representing the maximum
cross-correlation coefficients ($r_{\rm max}$). We used Monte Carlo
simulations incorporating flux randomization (FR) and random subset selection
(RSS) to calculate uncertainties in the time lags, thereby accounting for flux
errors and sampling uncertainties (Peterson et al., 1998). This resulted in a
cross-correlation centroid distribution (CCCD), where the 15.87% and 84.13%
quantiles represent the lower and upper bounds of time lags, respectively.
These time lags are tabulated in Table 2. Figure 1 presents the results of
auto-correlation functions (ACFs), CCFs, and CCCDs for a subset of targets
within the sample. It should be noted that for 5 objects, the emission line
lags cannot be reliably determined. Examination of the light curves reveals
that their emission line variability amplitudes are comparable to the errors
(Table 1), suggesting that it may be due to insufficient S/N in the spectra,
or potentially attributed to the rare anomalous response of the BLR, known as
an “AGN holiday” (Goad et al., 2016; Pei et al., 2017).

Furthermore, we conducted measurements of time delays for emission lines
across different velocity bins relative to continuum variations, referred to
as velocity-resolved delays. The process involved the following steps.
Initially, we generated the rms spectrum using the broad emission line
profiles obtained from the previous section. The emission line in the rms
spectrum was subsequently divided into several segments, each ensuring
reliable lag measurement by selecting smaller bins with sufficient variation
and containing identical integrated flux. Next, we applied these velocity bins
to each spectrum, yielding light curves for different velocity intervals. We
then separately calculated the time delays for H$\alpha$ and H$\beta$ within
each of these velocity bins, enabling an exploration of the relationship
between the BD and time lag across velocity space. The results are displayed
in Figure 3 and 4.

Figure 5: Velocity-resolved emission line flux ratios of Arp 151 and NGC 2617.
The black points represent the light curves of the continuum. The color maps
show the distributions of $F_{\rm H\alpha}/F_{\rm H\beta}$ values at different
velocities and times, where $F_{\rm H\alpha}/F_{\rm H\beta}$ values vary with
velocity and time. The black lines indicate the velocity-resolved time-
averaged values of $F_{\rm H\alpha}/F_{\rm H\beta}$, and the red dashed lines
represent the corresponding rms values.

###  3.3 Measurement of Emission Line Flux Ratios

Apart from determining time delays of emission lines, to understand the
relationship between the continuum and emission line ratios, we also computed
flux ratios between emission lines, such as $F_{\rm H\alpha}/F_{\rm H\beta}$
and $F_{\rm H\alpha}/F_{\rm H\gamma}$. To demonstrate the robustness of the
relationship, we calculated the flux ratios of the emission lines using three
distinct methods. The first results were derived from direct observations
carried out on the same night, with ratio errors estimated using the
uncertainties of spectral fitting. We used linear interpolation to align the
dates of these ratio measurements with the continuum light curve, resulting in
a newly adjusted continuum light curve. Considering the time delay between
H$\alpha$ and H$\beta$, this could weaken the correlation between the directly
measured BD and the continuum, as illustrated in Figure 1, where the minimum
value of the CCF correlation is delayed. To achieve a more accurate
correlation, correcting the lag of each emission line is necessary. To
mitigate the effects of both interpolation and time delay, we further examined
the correlation between flux ratios of emission lines and H$\beta$.

Direct observation results neither consider the influence of outliers in the
original light curves nor the difference in time delay. To test such factors
in our analysis, we utilized the code
JAVELIN333https://github.com/nye17/javelin, which models the continuum and
emission line variations using the DRW process (Zu et al., 2011), to
reconstruct the light curve. We derived emission line flux ratios and
uncertainties from these modeled curves, both at the same observation time and
after time delay correction. We then conducted a comparison of the line ratios
as determined by the direct observational method and JAVELIN-based methods,
which are depicted in Figure 1 (showing three AGNs examples). The relationship
remains robust across all three methods. Nevertheless, the presence of time
delays and outliers will contribute to a dispersion in the results, reducing
the apparent strength of the correlation. We also measured $\rm\Delta BD$,
which was calculated by taking the median of the five highest values and
subtracting the median of the five lowest values using the JAVELIN modeled
curves after time lag correction. These values are detailed in Table 1.

In addition, we measured the flux ratios of emission lines as a function of
velocity to investigate variation trend with line-of-sight velocity. We first
resampled broad H$\alpha$ and H$\beta$ profiles onto a shared velocity grid.
We then computed flux ratios from the velocity-binned line profiles at each
epoch. These velocity-dependent ratios are displayed against velocity and time
in Figure 5 to illustrate the more detailed variations. We also calculated
flux ratios of emission lines in velocity bins that correspond to the
velocity-resolved delays (as shown in the top panels of Figures 3 and 4). By
integrating the resampled H$\alpha$ and H$\beta$ profiles at each velocity
edge, we obtained light curves of $F_{\rm H\alpha}/F_{\rm H\beta}$ in these
bins. Finally, we derived the weighted average velocity-resolved flux ratios
by utilizing the $F_{\rm H\alpha}/F_{\rm H\beta}$ light curves at each
velocity edge. It is important to note that these ratios were not corrected
for the time delays of H$\alpha$ and H$\beta$ relative to the continuum. The
reason for this is that in our analysis, we used the average BD value for each
velocity bin, which can eliminate the influence of time delay (see tests in
Appendix A). To estimate the uncertainties of these ratios, we employed a
bootstrap approach, similar to the approach for CCF error calculations, with
10,000 iterations. This process accounted for errors in the emission line flux
ratios at each data point and the influence of the sampling. These results are
illustrated in Figures 3 and 4.

Figure 6: Simulated TFs for modeled BLR configurations. The upper-left panel
shows circular BLR with Keplerian motion. The upper-right panel shows circular
BLR with Keplerian and inflow motion. The lower-left panel shows circular BLR
with Keplerian and outflow motion. The lower-right panel shows elliptical BLR
with Keplerian motion. In each case, the red step line in the bottom panel
represents the corresponding steady-state emission line profile.

###  3.4 Transfer Function

According to the fundamental principles of the RM technique (e.g., Blandford &
McKee, 1982; Welsh & Horne, 1991; Bentz et al., 2010), the relationship
between the light curves of the continuum and an emission line can be
expressed by the following equation:

$L(\nu,t)=\int_{-\infty}^{\infty}\Psi(\nu,\tau)C(t-\tau)d\tau,$ (3)

where $\Psi(\nu,t)$ is the transfer function (TF) encoding the geometry and
kinematics of the BLR. The terms $L(\nu,t)$ and $C(t-\tau)$ stand for the
emission line and the continuum, respectively. Given the forms of the
continuum and the emission lines, we can determine the corresponding TF for a
specific BLR.

Studies on the BLRs of AGNs suggest that the BLR can be considered as a
collection of discrete clouds surrounding a single SMBH. The radial emissivity
distribution of these clouds can be well represented by a shifted
$\Gamma$-distribution (Pancoast et al., 2014):

$R_{\rm ga}=R_{\rm S}+FR_{\rm BLR}+g(1-F)\beta^{2}R_{\rm BLR},$ (4)

where $R_{\rm S}=2GM_{\rm{BH}}/c^{2}$ is the Schwarzschild radius of a black
hole with a mass $M_{\rm{BH}}$, $R_{\rm BLR}$ is the mean BLR radius, $F$ is
the ratio of the minimum radius to the mean BLR radius, $\beta$ is the shape
parameter, and g is randomly drawn from a $\Gamma$-distribution. For a
positive variable $x$, the $\Gamma$-distribution is described as:

$P(x|\alpha,\theta)\propto x^{\alpha-1}\rm exp(-x/\theta),$ (5)

where $\theta$ and $\alpha$ represent the scale parameter and shape parameter,
respectively. In this paper, we set $\beta$ =0.9, $F$ =0.6, $\theta$ = 1, and
$\alpha=1/\beta^{2}$ for simplicity. Additionally, we take into account the
angular distribution of the BLR clouds, which can be described by:

$\theta_{\rm mom}={\rm arccos}[{\rm cos}\,\theta_{\rm o}+(1-{\rm
cos}\,\theta_{\rm o})\times U^{\gamma}],$ (6)

where $\theta_{\rm mom}$ represents the angle between the angular momentum of
a cloud in the BLR relative to the BLR mid-plane. $U$ is a random number, and
$U^{\gamma}$ is considered to describe the concentration of BLR clouds
(Pancoast et al., 2014). In this work, we set $U^{\gamma}$ as a random number
between 0 and 1, indicating that the BLR clouds are uniformly distributed but
confined within $\pm\theta_{\rm o}$ range.

Here, we primarily model four disk-like BLRs with an opening angle of
$\theta_{\rm o}=$ 20 degrees. This configuration suggests that the clouds are
distributed within 20 degrees above and below the disk plane, representing a
flattened disk geometry that more accurately reflects the realistic BLR than
infinitely thin disk or spherical BLR geometry. Initially, we simulate the
Keplerian elliptical and circular BLR to test the effect of geometry, followed
by incorporating inflow and outflow motion into the Keplerian circular BLR.
For each case, simulations are conducted with 100,000 particles (clouds) under
a black hole mass of $2\times 10^{7}M_{\odot}$, similar to estimates for NGC
3516 (Feng et al., 2021a). The resulting TFs and steady-state emission line
profiles are presented in the following panels of Figure 6.

Table 3: Correlation measurements

| ${\rm con}~{}^{a}$ | ${\rm H\beta~{}^{b}}$
---|---|---
Object | $r_{\rm s}$ | $p_{\rm null}$ | $r_{\rm s}$ | $p_{\rm null}$
LAMP2008 campaign
Mrk 142 | 0.16 | 0.28 | -0.27 | 0.07
SBS 1116+583A | -0.34 | 0.03 | -0.74 | $2\times 10^{-08}$
Arp 151 | -0.84 | $5\times 10^{-11}$ | -0.80 | $5\times 10^{-10}$
Mrk 1310 | -0.54 | $2\times 10^{-04}$ | -0.83 | $6\times 10^{-13}$
Mrk 202 | -0.56 | $5\times 10^{-05}$ | -0.68 | $2\times 10^{-07}$
NGC 4253 | -0.21 | 0.19 | -0.77 | $2\times 10^{-10}$
NGC 4748 | -0.20 | 0.19 | -0.09 | 0.56
IC 4218 | -0.36 | 0.04 | -0.95 | $5\times 10^{-18}$
MCG -06-30-15 | -0.06 | 0.71 | -0.80 | $6\times 10^{-10}$
NGC 5548 | -0.49 | $6\times 10^{-04}$ | -0.68 | $1\times 10^{-07}$
Mrk 290 | -0.35 | 0.02 | -0.64 | $1\times 10^{-06}$
IC 1198 | -0.11 | 0.48 | -0.85 | $2\times 10^{-13}$
NGC 6814 | -0.56 | $1\times 10^{-04}$ | -0.63 | $5\times 10^{-06}$
LAMP2011 campaign
Mrk 40 (Arp 151) | -0.82 | $4\times 10^{-10}$ | -0.75 | $6\times 10^{-08}$
Mrk 50 | -0.93 | $4\times 10^{-19}$ | -0.94 | $8\times 10^{-20}$
Mrk 279 | -0.18 | 0.33 | -0.76 | $2\times 10^{-07}$
Mrk 704 | -0.13 | 0.45 | -0.09 | 0.59
Mrk 1511 | -0.63 | $2\times 10^{-05}$ | -0.62 | $4\times 10^{-05}$
NGC 4593 | -0.88 | $1\times 10^{-14}$ | -0.73 | $5\times 10^{-08}$
Zw 229-015 | -0.63 | $9\times 10^{-04}$ | -0.70 | $2\times 10^{-04}$
Lijiang CL-AGN campaign
NGC 2617 | -0.58 | $7\times 10^{-07}$ | -0.88 | $3\times 10^{-21}$
NGC 3516 | -0.67 | $7\times 10^{-10}$ | -0.92 | $9\times 10^{-28}$
NGC 4151 | -0.49 | $1\times 10^{-03}$ | -0.84 | $1\times 10^{-11}$

Note. — The objects analyzed in this study were selected from the LAMP2008,
LAMP2011, and Lijiang campaigns. We determined the Spearman rank correlation
between BD and continuum, as well as BD and H$\beta$. The results include the
Spearman rank correlation coefficient ($r_{\rm s}$) and probability of the
null hypothesis ($p_{\rm null}$).

a The Spearman rank correlation analysis between continuum and BD.

b The Spearman rank correlation analysis between H$\beta$ and BD.

##  4 Results

###  4.1 Anti-correlation between BD and Continuum

The BD is defined as the relative intensities of the hydrogen Balmer series
emission lines. H$\alpha$ and H$\beta$ are among the strongest optical
emission lines, and their flux ratio serves as the most frequently utilized
BD. We compiled a sample of 22 AGNs with simultaneous RM observations of both
H$\alpha$ and H$\beta$ from the LAMP2008 (Bentz et al., 2009) and LAMP2011
(Barth et al., 2015) and Lijiang CL-AGN RM campaign(Feng et al., 2021a, b; Li
et al., 2022), and some of these sources also included H$\gamma$ observations.
We conducted spectral fitting for 15 objects (from LAMP2008 and Lijiang
campaigns) with available spectral data to obtain pure broad emission line
profiles. For the remaining sources, we adopted the spectral decomposition
results from Barth et al. (2015). By utilizing the fluxes extracted from the
fitted broad emission lines, we can initiate an investigation into the
variability of the BD.

Figure 1 illustrates examples of BD light curves obtained from the LAMP2008,
LAMP2011, and Lijiang projects. It is evident that the BD can change by more
than 1 even within a timescale of $\sim$2 months, and similar variability
amplitudes are observed in most of the other sources. This suggests that the
BD in AGNs can exhibit considerable variability. Comparing the BD and
continuum light curves reveals almost completely opposite variability trends.
To quantify their correlation, we employed the Spearman rank correlation and
found a significant anti-correlation in 15 sources with 16 observations (see
Table 3). Therefore, decreasing BD with increasing continuum intensity appears
to be a common phenomenon. Upon further examination of the data, we noticed
that the remaining 7 sources, which did not show an inverse correlation, were
influenced by small variability amplitudes, non-synchronous observation times
between the continuum and emission line, and outliers. Furthermore, 5 of these
sources did not have a reliable time delay measurement from the CCF, because
their CCFs displayed anomalous profiles with $r_{\rm max}<$ 0.5. Notably, our
sample (including previous RM observation data) showed a significant time
delay between H$\alpha$ and H$\beta$ variations in many AGNs. Variability
introduces additional errors when attempting to measure the BD values at the
same incident photon level in the BLR using a single spectrum.

To assess the impact of these factors on our conclusions, we employed a method
based on the DRW model to fit multiple light curves simultaneously (see
details in Section 3.3). This approach can mitigate the influence of outliers
on individual light curves and enable us to test the effects of time delay on
correlation. We found a generally enhanced correlation between the continuum
and the BD measured with the DRW-based light curves, while almost all
theoretical DRW results corrected for lag display an inverse correlation.
Furthermore, by substituting H$\beta$ for the continuum, we can somewhat
reduce the effects of interpolation and lag, and found that 19 sources show a
significant inverse correlation. We also conducted tests using the BD values
measured from H$\alpha$ and H$\gamma$, and the results exhibit an identical
relationship. This relationship further supports the notion that the BD values
in the BLR are anti-correlated with the continuum fluxes.

The inverse relation between BD and the continuum can be explained by the
intrinsic Baldwin effect, which follows a power law relationship (e.g.,
Gilbert & Peterson, 2003), where $F_{\rm line}\propto F_{\rm con}^{m}$.
Specifically, for H$\alpha$ and H$\beta$, they can be expressed as powers of
$m1$ and $m2$, respectively. Then the BD is determined by $F_{\rm
H\alpha}/F_{\rm H\beta}\propto F_{\rm con}^{(m1-m2)}$, and will be anti-
correlated with the continuum when $m1<m2$, where $m$ represents the effective
line responsivity. This is consistent with the predictions from
photoionization models (Korista & Goad, 2004).

###  4.2 Results of Velocity-resolved IM and RM

The inherent ionization property leads to the natural inference that gas
clouds in the BLR located farther from the central ionizing source should
exhibit larger BD values. Therefore, we can study the gas distribution in the
radiation region by measuring the BD values at different velocities of broad
emission lines. In other words, comparing BD values between different velocity
bins allows us to obtain the relative distances of the emitting gas clouds
from the center. This technique is similar to velocity-resolved RM but is
independent of line-of-sight and kinematics. When combining these two methods
for analysis, we can effectively break the degeneracy of BLR geometry and
kinematics.

This kind of study requires high data quality, and we only analyzed 5 objects
that can simultaneously measure both H$\alpha$ and H$\beta$ velocity-resolved
lags. Initially, we examined the variation in the distribution of BD along the
velocity direction – the velocity-resolved BDs – using the fitted broad
emission line profiles. Figure 5 presents an example of BD as a function of
velocity and time for the LAMP2008 and Lijiang observations. As expected, we
do observe substantial variations in BD values at different velocities.
Notably, the BDs at any given velocity consistently exhibit an inverse
correlation with the continuum, and the overall shapes of the velocity-
resolved BDs are almost identical across each epoch. This inverse correlation
indicates that the BLR geometry should be constant throughout the monitoring
period. To mitigate the influence of variability, we constructed an average of
the BD values at each velocity across all observations. Such a process allowed
us to examine the ultimate form of the BD profile across different velocities.
Our findings revealed that the BD is symmetric around zero velocity in NGC
2617 and SBS 1116+583A. In contrast, Arp 151 and NGC 3516 display a decrease
in BD with increasing velocity, while NGC 4151 exhibits the opposite trend.

To compare the relationship between the BDs and lags at different line-of-
sight velocities of the broad emission line, we divided the rms spectrum of
broad H$\alpha$ into several bins of equal flux, and then measured the BD and
lag values within each bin. We found that, except for NGC 4151, the remaining
four sources exhibit a robust positive correlation between the BD and lag
values (see Figure 3 and 4). This correlation suggests that gas emitting with
longer time delays is indeed located at a considerable distance from the
ionizing source for these objects. The validity of this inference can be
further confirmed by analyzing the variation in the H$\alpha$-to-H$\beta$
ratio at different velocities. In most cases, where gas density exceeds
$10^{10}$ cm-3 and the $F_{\rm H\alpha}/F_{\rm H\beta}$ ratio falls within the
range of 2 to 10 (as shown in Figure 5 of Korista & Goad, 2004), a higher
hydrogen-ionizing photon flux (i.e., a smaller distance from the center)
corresponds to a less pronounced response in the BD. Specially, at
sufficiently high levels of ionizing photon flux, the BD becomes nearly
insensitive to continuum variations. Our approach primarily involved assessing
the rms values at different velocities to capture the BD variability, despite
the inherent uncertainties associated with rms due to its sensitivity to S/N
ratios. The distribution of the BD values is broadly consistent with the
distribution of their rms values, as shown in Figure 5. Intriguingly, Arp 151
and NGC 3516 exhibit distinct distances from the blue and red sides of the
emission line to the center, suggesting the existence of an asymmetric
geometry in the BLR. Moreover, the anti-correlation of BD and lag in NGC 4151
could serve as evidence of the complex kinematics within the BLR (Chen et al.,
2023).

##  5 Discussion

###  5.1 Mechanisms of BD Variation

It has been reported in some previous work that there is a negative
correlation between BD and continuum over years to decades in several
individual AGNs (Shapovalova et al., 2004; Popović et al., 2011; Ma et al.,
2023; Oknyansky et al., 2023). This phenomenon, especially when accounting for
the CL mechanisms, has been partly ascribed to dust extinction variations
(Kollatschny et al., 2000), whilst others propose ionization effects within
the BLR as the principal cause (Shapovalova et al., 2010). It is generally
believed that dust grains are located beyond the sublimation radius (around
3-5 times the BLR radius), otherwise they would be vaporized within hours
(Baskin & Laor, 2018). Changes in dust clouds on these size scales typically
require decades (LaMassa et al., 2015), which could explain some long-term BD
variations. However, the objects in our selected sample were monitored for
only 2-6 months, and the variability timescales of their BDs are even shorter
at $\sim$ 2 months. This discrepancy suggests that it is implausible to expect
significant changes in the obscuring material. To validate our hypothesis, we
conducted two tests (see Appendix B for details): First, we used Zwicky
Transient Facility (ZTF) data to examine variability of spectral index for
these AGNs (often used as an extinction proxy), and found that their spectral
indices remained nearly constant on a timescale of $\sim$5 years. Second, we
simulated the BD-continuum relationship under varying extinction conditions,
which did not reproduce the observed patterns, further excluding the impact of
extinction.

The detection of delayed emission line responses relative to continuum changes
for all substantially variable sources contradicts expectations from
extinction models (see Table 2). Unified AGN models posit a dusty torus
surrounding the outer BLR (Antonucci, 1993). If extinction changes drove
variations, the BLR flux would respond first before the accretion disk
continuum. This would result in emission line variations preceding continuum
changes, which are not observed. Therefore, the anti-correlation between BD
and the continuum shown in our study reflects an ionization effect rather than
extinction. For a finite BLR size, the smaller mean formation radius and
larger effective responsivity of H$\beta$ will result in a lower BD at
increased flux levels. In reality, an increase in continuum emission can
illuminate more distant BLR regions, which enhances the broad Balmer emission
lines with larger local BDs. Despite this, the overall BD measured across the
BLR is observed to decline, primarily due to the overwhelming contribution
from closely orbiting, more heavily-ionized inner clouds, which possess lower
intrinsic BDs. This effect becomes more pronounced with increasing central
ionizing flux. The observed BD-continuum anti-correlation reinforces the
fundamental relationship between local BD and radius, which is necessary to
the velocity-resolved IM technique.

Figure 7: A cartoon depiction of kinematics and geometry effects on velocity-
resolved BDs and lags. The upper-left panel shows circular BLR with Keplerian
motion. The upper-right shows circular BLR with Keplerian and inflow motion.
The lower-left shows circular BLR with Keplerian and outflow motion. The
lower-right shows elliptical BLR with Keplerian motion. In each case, the
dashed line is the trajectory of a cloud, and the arrows are the motion
direction of the cloud. The green line is the direction of the photon
propagation. The red color represents the part of BLR where emission lines are
redward shifted with respect to the black hole, and the blue color represents
the part of BLR where emission lines are blueshifted.

###  5.2 Velocity-Resolved IM Technique

Velocity-resolved IM technique assumes that these local BDs of BLR gas clouds
increase monotonically with distance from the ionization source, a reasonable
premise for idealized dust-free environments with homogeneous matter
distributions. Real AGNs deviate from these idealized conditions - the
ubiquitous presence of dust can impact observed BDs through internal
extinction (Czerny & Hryniewicz, 2011; Baskin & Laor, 2018). Radial dust
density gradients necessitate re-evaluating the local BD-distance correlation
even if the overall extinction level holds steady. Near the outer boundary of
BLR where it transitions into the dusty torus, mixing of gas and dust becomes
increasingly likely as temperatures approach dust sublimation thresholds.
Here, outwardly rising dust densities could amplify BDs beyond these
increments predicted by dust-free models.

The ionization state of the BLR is not solely determined by the incident
continuum, but is also intricately linked to the physical condition of the gas
(Ilić et al., 2012). When significant changes occur in the physical parameters
within the BLR, the BD value will also vary accordingly. One key parameter
influencing the positive BD-distance relationship is hydrogen density. In the
vicinity of the central region, where gas density is high, this leads to an
increased tendency for collisional excitation, which enhances the emission of
the H$\alpha$ line and thereby increases the BD value. Meanwhile, as gas
density reaches a certain threshold, the BD value tends to stabilize or even
inversely correlate with density increases (Korista & Goad, 2004). This is
attributed to the reduced efficiency of converting ionizing photons into
escaping line photons at high gas densities, particularly affecting H$\alpha$
with its higher optical depth compared to H$\beta$. These properties may lead
to a more complex distribution of BD values within the BLR, where the former
effect seems to contradict the foundational premises of velocity-resolved IM
and introduces some uncertainty in its application, while the latter effect
still supports the assumption of BD increasing with distance.

Current studies indicate a weak correlation between BD and hydrogen density.
Even significant changes in density can induce only modest changes in BD, as
illustrated by the fact that a density contrast of several orders of magnitude
between the BLR and the narrow line region typically translates into a
predicted BD variation of about 0.25 under Case B recombination conditions.
Our results demonstrate that a 20% fluctuation in the continuum can result in
a significant change in the BD, with absolute values greater than 1 in some
objects as shown in Figure 1 and Table 1. This indicates that within the BLR,
where hydrogen density is expected to decrease with distance following a
power-law of $n(R)\propto R^{-s}$ ($s$ often between 0 and 2)—the effect of
the incident continuum on BD is more pronounced than the effect of varying
hydrogen densities. Given that the incident continuum flux diminishes
approximately as $F(R)\propto R^{-2}$, outpacing the decrease in hydrogen
density, we conclude that the ionizing effect of the continuum is the dominant
factor in determining the BD in AGNs. In reality, the presence of dust or
higher gas density in certain regions of the BLR (Netzer, 2015) can lead to a
more rapid decline in the incident continuum flux $F(R)$ due to absorption and
scattering (Baskin & Laor, 2018), further reinforcing our conclusion.

Photoionization simulations show that the emissivity distributions of
H$\alpha$ and H$\beta$ in BLR do exhibit differences (or different transfer
functions), displaying a trend that $F_{\rm H\beta}$ decreases faster than
$F_{\rm H\alpha}$ with the increasing distance from the center (Guo et al.,
2020; Wu et al., 2023). The implications of this property are significant for
velocity-resolved IM, as it suggests the BD value is a monotonic function of
radius. For a symmetric Keplerian BLR (Horne et al., 2004), the detected BD
and lag values should be consistent for the emitting gas clouds with the same
line-of-sight velocities on the blue and red sides of broad emission line
because these clouds are equidistant to the center of the BLR. Considering
that the gas cloud closer to the central black hole moves at a higher velocity
and receives a stronger flux of the incident continuum, the detected lag is
expected to be shorter, and the BD value should be smaller (see discussion in
Korista & Goad, 2004, for more details). A decreasing/increasing trend of lag
from the blue side to the red side of broad emission line is a key
characteristic previously attributed to inflow/outflow within the BLR. In this
scenario, the BLR is still considered to be symmetric, but the viewing angle
causes the blue/red-shifted gas to appear more distant, resulting in a longer
lag.

However, for a symmetric BLR, regardless of its dynamical state, the BD
profile should maintain symmetry in the direction of line-of-sight velocity.
This is due to the fact that the BD value is solely dependent on its distance
from the central ionizing source and is independent of the propagation path of
the light. Figure 7 gives a schematic cartoon showing the effects of
kinematics on the velocity-resolved lags and BDs, and the corresponding
simulated two-dimensional velocity-resolved time lag maps are shown in Figure
6. These simulations illustrate that both inflow and outflow can result in
asymmetric velocity-resolved lag profiles across the broad emission line,
assuming that the ionization properties are not influenced by kinematics.

###  5.3 Diverse geometry and kinematics of BLR

Our results reveal remarkable diversity in the geometry and kinematics of the
BLR among different AGNs. For NGC 2617 and SBS 1116+583A, the symmetric nature
of the velocity-resolved BD profiles, together with corresponding time lags,
are indicative of a BLR that predominantly exhibits Keplerian motion in a
symmetric geometry (Figure 7). Conversely, the scenarios presented by Arp 151,
NGC 3516, and NGC 4151 imply a starkly different picture. Their longer lags in
the blue wing suggest the possibility of inflowing characteristics or the
blue-shifted gas being located farther from the center. In the context of Arp
151 and NGC 3516, distinct BD-lag correlations on the blue and red wings point
to non-axisymmetric geometry.

Intriguingly, NGC 4151 shows an anti-correlation between velocity-resolved BD
and lag values. This peculiar signature may be attributed to special motions
within the BLR. Indeed, we noticed that the shape of the velocity-resolved
H$\beta$ lags exhibited temporal changes during the 2020-2021 period (Chen et
al., 2023), but the results from additional epochs between 2018 and 2022 align
with our velocity-resolved BDs conducted in 2020-2021. The transient kinematic
shifts could be related to the dramatic fluctuations in radiation pressure
owing to changes in AGN luminosity, or some localized perturbations within the
BLR, neither of which should substantially impact on the global geometry of
the BLR. It must be emphasized that the aforementioned conclusions are
predicated on the assumption of symmetrical kinematics within the BLR. The
presence of an asymmetric inflow or outflow could also replicate similar
observational signatures under finely-tuned viewing angles. For instance, in
Arp 151 and NGC 3516, we might coincidentally detect only the blueshifted
inflowing gas, whereas, in NGC 4151, the situation could be the inverse. These
kinematic properties would likely be inconsistent with the virialized motions
and can be tested through multiple RM observations.

There are several potential configurations of geometry for such an asymmetric
BLR. One possible scenario is that the BLR contains an uneven spiral arm or
hot spot, creating a stronger response at some particular velocity (Newman et
al., 1997; Storchi-Bergmann et al., 2003). These substructures can change as
the BLR undergoes rapid motion, causing variations in the velocity-resolved
signal as shown in NGC 4151. The velocity-resolved RM measurements of Arp 151
and NGC 3516 have remained virtually unchanged for several years (Denney et
al., 2010; Pancoast et al., 2018; Feng et al., 2021a), requiring a BLR with
stable geometry and kinematics. A Keplerian elliptical disk model may satisfy
this condition (Eracleous et al., 1995). This configuration can persist for
long durations and will result in unequal distances between the red and blue
sides of the BLR to the central ionizing source (Figure 7). Also, gas clouds
in elliptical orbits would have orbital velocities varying with their
positions. This would lead to an observational asymmetric profile of broad
emission line due to the longer orbiting time at the apocenter. Figure 6 shows
the transfer functions and emission line profiles given by our simulations of
elliptical disk BLR. It is evident that the emission line on the blue side is
noticeably stronger than on the red side, and the lag on the blue side is
generally larger than on the red side. These findings are consistent with
observational results in NGC 3516. Binary SMBH can also lead to asymmetric BLR
(Eracleous et al., 1997), but this situation is generally believed to produce
periodic variability, which are rarely detected reliably in AGNs.

###  5.4 Implications for AGN physics

The architecture of an asymmetric BLR might be related to the formation and
evolution of AGNs, such as the merging of SMBHs, the formation of BLR, the
supply of matter to the accretion disk, and the feedback mechanisms. Eracleous
et al. (1995) proposed that an elliptical disk in the BLR could be associated
with binary black holes or the tidal disruption of stars by an SMBH. Wang et
al. (2017) suggested that tidally disrupted clumps by the central SMBH from
the dusty torus can form the BLR, including spiral-in gas as inﬂow,
circularized gas, and ejecta as outﬂow. In this framework, the evolutionary
timescale of the BLR typically requires decades or longer, consistent with the
observed stability of the BLR in Arp 151, NGC 3516. The detection of a tidal
disruption event in 1ES 1927+654 also supports this theory, which suggests
that the BLR originates from the torus (Trakhtenbrot et al., 2019). If the
inward spiraling process is disturbed, for instance, by a smaller black hole
or a star passing through, it could form a spiral arm structure or an
elliptical disk. Alternatively, the BLR could simply be a result of failed
radiation-driven dusty outflows from the accretion disk (Czerny et al., 2017).
In this model, if the powerful winds are anisotropic, it can also lead to an
asymmetric BLR. Furthermore, when there are enough clouds in the BLR, their
self-gravity can also produce some non-uniform structures (Du & Wang, 2023).

While previous observations have suggested the possibility of asymmetry within
the BLR, and various theoretical models have been proposed to explain this
asymmetric geometry, there is still a tendency to utilize simpler symmetrical
models when interpreting observational results. The joint analysis of
velocity-resolved BDs and time lags can effectively reveal non-axisymmetric
geometries of the BLR. In the future, as we confirm more instances of such
asymmetry, it will enable us to investigate the correlation between the
geometry of the BLR and the physical parameters of AGN, such as the accretion
rate and the black hole mass. These relationships are vital to our
comprehension of the formation of the BLR and the internal physical processes
in AGNs. Furthermore, the geometry of BLR is essential to accurately measure
$M_{\rm{BH}}$, as emission line profiles and time delays resulting from
asymmetric BLRs significantly differ from expectations of circular disk-shaped
BLR models (Figure 6). Therefore, further research into asymmetric BLRs is
warranted to improve our understanding of AGN dynamics and measurements of
$M_{\rm{BH}}$.

##  6 Summary

In this work, we utilized high-quality spectroscopic monitoring data from the
LAMP2008, LAMP2011, and Lijiang projects to conduct a detailed analysis of the
geometric and kinematic characteristics of the BLR in AGNs. Specifically, we
examined the flux ratios and time delays of the H$\alpha$ and H$\beta$
emission lines in this sample. By integrating the new technique of velocity-
resolved IM with traditional velocity-resolved RM, we unveiled the complex
kinematics and geometry within the BLR. Our key results are as follows.

  1. 1.

We detected a significant anti-correlation between the BD and continuum on
timescales of weeks to months in 15 out of 22 AGNs. Dust extinction variations
were excluded as the primary driver, indicating that this phenomenon arises
from ionization effects within the BLR. The decreasing BD with increasing
ionizing flux is consistent with photoionization model predictions.

  2. 2.

Tight correlations between velocity-resolved BDs and RM lags were found in 4
out of 5 AGNs with high-quality data. The BD profiles of NGC 2617 and SBS
1116+583A are symmetric around zero velocity, suggesting a Keplerian BLR with
symmetric geometry. In contrast, the asymmetric BD profiles and distinct lags
on the blue and red sides in Arp 151 and NGC 3516 point to non-axisymmetric
BLR geometry, possibly an elliptical disk or spiral arm structure.

  3. 3.

NGC 4151 exhibits an anti-correlation between velocity-resolved BDs and lags,
implying complex kinematics within the BLR. However, its BD profile remains
stable over 4 years despite transient kinematic changes, supporting a
persistent asymmetric geometry.

  4. 4.

We discussed the assumptions and limitations of the IM technique. While
factors like radial dust distribution and gas density may affect local BDs,
the incident continuum flux is likely the dominant driver of the BD-distance
relation within the BLR. Simulations of various BLR configurations
demonstrated the power of combining IM and RM to distinguish between geometric
and kinematic effects.

  5. 5.

The detection of asymmetric BLR geometry has important implications for AGN
physics. It may represent a formation imprint related to SMBH mergers, tidally
disrupted clumps from the torus, or radiation-driven disk winds. Asymmetry can
also impact SMBH mass measurements. Further investigations into the
correlations between BLR geometry and AGN properties are warranted.

This study showcases the potential of the velocity-resolved IM technique in
probing the BLR geometry and its promise in advancing our understanding of AGN
physics when combined with RM. Future applications to larger samples will
provide valuable insights into the formation and evolution of AGNs and refine
SMBH mass estimates. As spectroscopic monitoring data continues to accumulate,
the IM method can be extended to more objects, opening up new avenues for
exploring the diversity of BLR geometry and kinematics across the AGN
population.

We appreciate the referee for constructive comments and valuable suggestions,
that greatly improved this paper. We are grateful to J.-R. Mao, J.-C. Wu, and
H.-X. Guo for the comments that improved the manuscript. This work is
supported by National Key R&D Program of China (No. 2021YFA1600404), the
National Natural Science Foundation of China (grants No. 12303022, 12203096,
12373018, and 11991051), Yunnan Fundamental Research Projects (grants NO.
202301AT070358 and 202301AT070339), Yunnan Postdoctoral Research Foundation
Funding Project, Special Research Assistant Funding Project of Chinese Academy
of Sciences, and the science research grants from the China Manned Space
Project with No. CMS-CSST-2021-A06.

We acknowledge the support of the staff of the Lijiang 2.4 m telescope.
Funding for the telescope has been provided by Chinese Academy of Sciences and
the People’s Government of Yunnan Province. This work has made use of data
from the Lick AGN Monitoring Project public data release. The Zwicky Transient
Facility Collaboration is supported by U.S. National Science Foundation
through the Mid-Scale Innovations Program (MSIP). This research has made use
of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet
Propulsion Laboratory, California Institute of Technology, under contract with
NASA.

##  Appendix A Examining the Impact of Reverberation Effect on Velocity-
resolved BD

To examine the reverberation effect in velocity-resolved flux ratios (BDs), we
conducted the relevant tests. Initially, we measured the light curves of
H$\alpha$ and H$\beta$ within individual velocity bins. The direct flux ratios
were determined from these light curves, and the average BDs for each bin are
depicted as blue points in Figure 8. Subsequently, we analyzed the time lags
of H$\alpha$ and H$\beta$ in each velocity bin relative to the continuum using
the CCF method. This allowed us to correct the time lags of each emission line
relative to the continuum in each velocity bin, resulting in lag-corrected
emission line light curves. The corresponding flux ratios were then calculated
in the same manner as before, depicted as orange points in Figure 8. Our
findings revealed that the average BD results and trends remained largely
consistent between the two methods.

Figure 8: Velocity-resolved emission line flux ratios for Arp 151, SBS
1116+583A, NGC 2617, NGC 3516 and NGC 4151. In each panel, the blue points
represent the flux ratios directly measured from H$\alpha$ and H$\beta$ light
curves within each velocity bin. The orange points show the flux ratios of
light curves that have been adjusted for the time lags of H$\alpha$ and
H$\beta$ relative to the continuum within each velocity bin. The vertical
dashed lines indicate the boundaries of the velocity bins, same as shown in
Figures 3 and 4.

##  Appendix B Examination of Dust Extinction Variability

Figure 9: The flux-to-flux diagrams between ZTF-g and ZTF-r bands. The data
samples are from LAMP2008, LAMP2011, and Lijiang projects. 19 out of 22 AGNs
were included in the analysis (3 AGNs lacking data). The colored points
represent different observational times, and the grey line represents the
linear fitting result.

Here, we implemented two distinct approaches to investigate potential changes
in the reddening within AGNs. The first approach utilizes flux-flux diagrams
constructed from multi-band light curves. These diagrams enable us to measure
the slopes between light curves in two bands, providing a means to assess
extinction as outlined by Winkler (1997). By tracking the evolution of these
slopes over time, we are able to identify any shifts in reddening assuming
that the AGN power-law index is constant.

For this initial approach, we utilized observational data from the ZTF, which
conducts time-domain observations using a 48-inch telescope across the ZTF-g,
ZTF-r, and ZTF-i bands (Bellm et al., 2019). Due to the relatively sparse data
in the ZTF-i band, we focused our analysis on light curves from the ZTF-g and
ZTF-r bands. We paired ZTF-g and ZTF-r band light curves with time intervals
of less than 0.8 days and converted them to fluxes based on the AB magnitude
system. Then, we corrected for Galactic extinction based on dust maps
(Schlafly & Finkbeiner, 2011) and band-specific extinction coefficients (Yang
et al., 2023). Using the Least Squares Method, we fitted these paired light
curves with a linear regression. As shown in Figure 9, the highly linear
relationship and lack of curvature in the flux-flux diagram for 19 AGNs (with
3 lacking sufficient data) suggests negligible reddening variations over
month-long timescales.

The second method relies on data from LAMP2008 and Lijiang samples, with the
LAMP2011 dataset contributing narrow line components to the light curves of
H$\alpha$ and H$\beta$. Under the assumption that continuum variations are
driven by fluctuating dust extinction, we utilized the extinction curve from
Fitzpatrick (1999) to estimate the corresponding reddening and to model the
expected $F_{\rm H\alpha}/F_{\rm H\beta}$. The results are illustrated in
Figure 10, where the red lines represent our simulated expectations. The poor
agreement between modeled and observed $F_{\rm H\alpha}/F_{\rm H\beta}$
(Figure 10) implies that variations in extinction alone may not explain the
observed changes in the continuum and the variability of $F_{\rm
H\alpha}/F_{\rm H\beta}$.

Figure 10: The simulated result between emission line flux ratio and
continuum. The data employed for this analysis originate from LAMP2008 and
Lijiang projects. In the figures, the observed data are depicted as black
points. Operating under the assumption that changes in the continuum result
from dust extinction fluctuations, we have modeled the expected $F_{\rm
H\alpha}/F_{\rm H\beta}$, which is represented by the red lines.

## References

  * Antonucci (1993) Antonucci, R. 1993, ARA&A, 31, 473. doi:10.1146/annurev.aa.31.090193.002353
  * Bao et al. (2022) Bao, D.-W., Brotherton, M. S., Du, P., et al. 2022, ApJS, 262, 14. doi:10.3847/1538-4365/ac7beb
  * Barth et al. (2015) Barth, A. J., Bennert, V. N., Canalizo, G., et al. 2015, ApJS, 217, 26. doi:10.1088/0067-0049/217/2/26
  * Baskin & Laor (2018) Baskin, A. & Laor, A. 2018, MNRAS, 474, 1970. doi:10.1093/mnras/stx2850
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002. doi:10.1088/1538-3873/aaecbe
  * Bentz et al. (2009) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2009, ApJ, 705, 199. doi:10.1088/0004-637X/705/1/19
  * Bentz et al. (2010) Bentz, M. C., Walsh, J. L., Barth, A. J., et al. 2010, ApJ, 716, 993. doi:10.1088/0004-637X/716/2/993
  * Blandford & McKee (1982) Blandford, R. D. & McKee, C. F. 1982, ApJ, 255, 419. doi:10.1086/159843
  * Boroson & Green (1992) Boroson, T. A. & Green, R. F. 1992, ApJS, 80, 109. doi:10.1086/191661
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000. doi:10.1046/j.1365-8711.2003.06897.x
  * Cackett & Horne (2006) Cackett, E. M. & Horne, K. 2006, MNRAS, 365, 1180. doi:10.1111/j.1365-2966.2005.09795.x
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245. doi:10.1086/167900
  * Chen & Halpern (1989) Chen, K. & Halpern, J. P. 1989, ApJ, 344, 115. doi:10.1086/167782
  * Chen et al. (1989) Chen, K., Halpern, J. P., & Filippenko, A. V. 1989, ApJ, 339, 742. doi:10.1086/167332
  * Chen et al. (2023) Chen, Y.-J., Bao, D.-W., Zhai, S., et al. 2023, MNRAS, 520, 1807. doi:10.1093/mnras/stad051
  * Clavel et al. (1991) Clavel, J., Reichert, G. A., Alloin, D., et al. 1991, ApJ, 366, 64. doi:10.1086/169540
  * Czerny & Hryniewicz (2011) Czerny, B. & Hryniewicz, K. 2011, A&A, 525, L8. doi:10.1051/0004-6361/201016025
  * Czerny et al. (2017) Czerny, B., Li, Y.-R., Hryniewicz, K., et al. 2017, ApJ, 846, 154. doi:10.3847/1538-4357/aa8810
  * Denney et al. (2010) Denney, K. D., Peterson, B. M., Pogge, R. W., et al. 2010, ApJ, 721, 715. doi:10.1088/0004-637X/721/1/715
  * Dong et al. (2008) Dong, X., Wang, T., Wang, J., et al. 2008, MNRAS, 383, 581. doi:10.1111/j.1365-2966.2007.12560.x
  * Du & Wang (2023) Du, P. & Wang, J.-M. 2023, A&A, 671, A26. doi:10.1051/0004-6361/202244780
  * Edelson et al. (2002) Edelson, R., Turner, T. J., Pounds, K., et al. 2002, ApJ, 568, 610. doi:10.1086/323779
  * Elitzur et al. (2014) Elitzur, M., Ho, L. C., & Trump, J. R. 2014, MNRAS, 438, 3340. doi:10.1093/mnras/stt2445
  * Emmering et al. (1992) Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, ApJ, 385, 460. doi:10.1086/170955
  * Eracleous et al. (1997) Eracleous, M., Halpern, J. P., M. Gilbert, A., et al. 1997, ApJ, 490, 216. doi:10.1086/304859
  * Eracleous et al. (1995) Eracleous, M., Livio, M., Halpern, J. P., et al. 1995, ApJ, 438, 610. doi:10.1086/175104
  * Feng et al. (2021a) Feng, H.-C., Hu, C., Li, S.-S., et al. 2021a, ApJ, 909, 18. doi:10.3847/1538-4357/abd851
  * Feng et al. (2021b) Feng, H.-C., Liu, H. T., Bai, J. M., et al. 2021b, ApJ, 912, 92. doi:10.3847/1538-4357/abefe0
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63. doi:10.1086/316293
  * Flohic et al. (2012) Flohic, H. M. L. G., Eracleous, M., & Bogdanović, T. 2012, ApJ, 753, 133. doi:10.1088/0004-637X/753/2/133
  * Gaskell & Peterson (1987) Gaskell, C. M. & Peterson, B. M. 1987, ApJS, 65, 1. doi:10.1086/191216
  * Gaskell & Sparke (1986) Gaskell, C. M. & Sparke, L. S. 1986, ApJ, 305, 175. doi:10.1086/164238
  * Gilbert et al. (1999) Gilbert, A. M., Eracleous, M., Filippenko, A. V., et al. 1999, Structure and Kinematics of Quasar Broad Line Regions, 175, 189
  * Gilbert & Peterson (2003) Gilbert, K. M. & Peterson, B. M. 2003, ApJ, 587, 123. doi:10.1086/368112
  * Goad et al. (2016) Goad, M. R., Korista, K. T., De Rosa, G., et al. 2016, ApJ, 824, 11. doi:10.3847/0004-637X/824/1/1
  * Goad et al. (2004) Goad, M. R., Korista, K. T., & Knigge, C. 2004, MNRAS, 352, 277. doi:10.1111/j.1365-2966.2004.07919.x
  * Goad et al. (2012) Goad, M. R., Korista, K. T., & Ruff, A. J. 2012, MNRAS, 426, 3086. doi:10.1111/j.1365-2966.2012.21808.x 1
  * Goad & Wanders (1996) Goad, M. & Wanders, I. 1996, ApJ, 469, 113. doi:10.1086/177764
  * Gravity Collaboration et al. (2018) Gravity Collaboration, Sturm, E., Dexter, J., et al. 2018, Nature, 563, 657. doi:10.1038/s41586-018-0731-9
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Amorim, A., Bauböck, M., et al. 2020, A&A, 643, A154. doi:10.1051/0004-6361/202039067
  * GRAVITY Collaboration et al. (2021) GRAVITY Collaboration, Amorim, A., Bauböck, M., et al. 2021, A&A, 648, A117. doi:10.1051/0004-6361/202040061
  * Grier et al. (2013) Grier, C. J., Peterson, B. M., Horne, K., et al. 2013, ApJ, 764, 47. doi:10.1088/0004-637X/764/1/47
  * Guo et al. (2020) Guo, H., Shen, Y., He, Z., et al. 2020, ApJ, 888, 58. doi:10.3847/1538-4357/ab5db0
  * Horne et al. (2004) Horne, K., Peterson, B. M., Collier, S. J., et al. 2004, PASP, 116, 465. doi:10.1086/420755
  * Hu et al. (2016) Hu, C., Wang, J.-M., Ho, L. C., et al. 2016, ApJ, 832, 197. doi:10.3847/0004-637X/832/2/197
  * Ilić et al. (2012) Ilić, D., Popović, L. Č., La Mura, G., et al. 2012, A&A, 543, A142. doi:10.1051/0004-6361/201219299
  * Kelly et al. (2009) Kelly, B. C., Bechtold, J., & Siemiginowska, A. 2009, ApJ, 698, 895. doi:10.1088/0004-637X/698/1/895
  * Kollatschny et al. (2000) Kollatschny, W., Bischoff, K., & Dietrich, M. 2000, A&A, 361, 901
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511. doi:10.1146/annurev-astro-082708-101811
  * Korista & Goad (2004) Korista, K. T. & Goad, M. R. 2004, ApJ, 606, 749. doi:10.1086/383193
  * LaMassa et al. (2015) LaMassa, S. M., Cales, S., Moran, E. C., et al. 2015, ApJ, 800, 144. doi:10.1088/0004-637X/800/2/144
  * Laor et al. (2006) Laor, A., Barth, A. J., Ho, L. C., et al. 2006, ApJ, 636, 83. doi:10.1086/497908
  * Li et al. (2022) Li, S.-S., Feng, H.-C., Liu, H. T., et al. 2022, ApJ, 936, 75. doi:10.3847/1538-4357/ac8745
  * Li et al. (2021) Li, S.-S., Yang, S., Yang, Z.-X., et al. 2021, ApJ, 920, 9. doi:10.3847/1538-4357/ac116e
  * Li et al. (2014) Li, Y.-R., Wang, J.-M., Hu, C., et al. 2014, ApJ, 786, L6. doi:10.1088/2041-8205/786/1/L6
  * Lu et al. (2022) Lu, K.-X., Bai, J.-M., Wang, J.-M., et al. 2022, ApJS, 263, 10. doi:10.3847/1538-4365/ac94d3
  * Lu et al. (2019) Lu, K.-X., Bai, J.-M., Zhang, Z.-X., et al. 2019, ApJ, 887, 135. doi:10.3847/1538-4357/ab5790
  * Ma et al. (2023) Ma, Y.-S., Li, S.-J., Gu, C.-S., et al. 2023, MNRAS, 522, 5680. doi:10.1093/mnras/stad1377
  * Mejía-Restrepo et al. (2018) Mejía-Restrepo, J. E., Lira, P., Netzer, H., et al. 2018, Nature Astronomy, 2, 63. doi:10.1038/s41550-017-0305-z
  * Netzer (2015) Netzer, H. 2015, ARA&A, 53, 365. doi:10.1146/annurev-astro-082214-122302
  * Newman et al. (1997) Newman, J. A., Eracleous, M., Filippenko, A. V., et al. 1997, ApJ, 485, 570. doi:10.1086/304460
  * Nicastro (2000) Nicastro, F. 2000, ApJ, 530, L65. doi:10.1086/312491
  * O’Donnell (1994) O’Donnell, J. E. 1994, ApJ, 422, 158. doi:10.1086/173713
  * Oknyansky et al. (2023) Oknyansky, V. L., Brotherton, M. S., Tsygankov, S. S., et al. 2023, MNRAS, 525, 2571. doi:10.1093/mnras/stad2442
  * Pancoast et al. (2018) Pancoast, A., Barth, A. J., Horne, K., et al. 2018, ApJ, 856, 108. doi:10.3847/1538-4357/aab3c6
  * Pancoast et al. (2014) Pancoast, A., Brewer, B. J., Treu, T., et al. 2014, MNRAS, 445, 3073. doi:10.1093/mnras/stu1419
  * Pei et al. (2017) Pei, L., Fausnaugh, M. M., Barth, A. J., et al. 2017, ApJ, 837, 131. doi:10.3847/1538-4357/aa5eb1
  * Peterson et al. (1998) Peterson, B. M., Wanders, I., Horne, K., et al. 1998, PASP, 110, 660. doi:10.1086/316177
  * Popović et al. (2011) Popović, L. Č., Shapovalova, A. I., Ilić, D., et al. 2011, A&A, 528, A130. doi:10.1051/0004-6361/201016317
  * Reynolds et al. (1997) Reynolds, C. S., Ward, M. J., Fabian, A. C., et al. 1997, MNRAS, 291, 403. doi:10.1093/mnras/291.3.403
  * Rodríguez-Pascual et al. (1997) Rodríguez-Pascual, P. M., Alloin, D., Clavel, J., et al. 1997, ApJS, 110, 9. doi:10.1086/312996
  * Salpeter (1964) Salpeter, E. E. 1964, ApJ, 140, 796. doi:10.1086/147973
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103. doi:10.1088/0004-637X/737/2/103
  * Shapovalova et al. (2004) Shapovalova, A. I., Doroshenko, V. T., Bochkarev, N. G., et al. 2004, A&A, 422, 925. doi:10.1051/0004-6361:20035652
  * Shapovalova et al. (2010) Shapovalova, A. I., Popović, L. Č., Burenkov, A. N., et al. 2010, A&A, 509, A106. doi:10.1051/0004-6361/200912311
  * Storchi-Bergmann et al. (2003) Storchi-Bergmann, T., Nemmen da Silva, R., Eracleous, M., et al. 2003, ApJ, 598, 956. doi:10.1086/378938
  * Trakhtenbrot et al. (2019) Trakhtenbrot, B., Arcavi, I., MacLeod, C. L., et al. 2019, ApJ, 883, 94. doi:10.3847/1538-4357/ab39e4
  * van Groningen & Wanders (1992) van Groningen, E. & Wanders, I. 1992, PASP, 104, 700. doi:10.1086/133039
  * Villafaña et al. (2022) Villafaña, L., Williams, P. R., Treu, T., et al. 2022, ApJ, 930, 52. doi:10.3847/1538-4357/ac6171
  * Walsh et al. (2009) Walsh, J. L., Minezaki, T., Bentz, M. C., et al. 2009, ApJS, 185, 156. doi:10.1088/0067-0049/185/1/156
  * Wanders et al. (1995) Wanders, I., Goad, M. R., Korista, K. T., et al. 1995, ApJ, 453, L87. doi:10.1086/309750
  * Wang et al. (2017) Wang, J.-M., Du, P., Brotherton, M. S., et al. 2017, Nature Astronomy, 1, 775. doi:10.1038/s41550-017-0264-4
  * Wang et al. (2020) Wang, S., Shen, Y., Jiang, L., et al. 2020, ApJ, 903, 51. doi:10.3847/1538-4357/abb36d
  * Waters et al. (2016) Waters, T., Kashi, A., Proga, D., et al. 2016, ApJ, 827, 53. doi:10.3847/0004-637X/827/1/53
  * Welsh & Horne (1991) Welsh, W. F. & Horne, K. 1991, ApJ, 379, 586. doi:10.1086/170530
  * Winkler (1997) Winkler, H. 1997, MNRAS, 292, 273. doi:10.1093/mnras/292.2.273
  * Wu et al. (2023) Wu, J., Wu, Q., Xue, H., et al. 2023, ApJ, 950, 106. doi:10.3847/1538-4357/acce9e
  * Yang et al. (2023) Yang, Q., Green, P. J., MacLeod, C. L., et al. 2023, ApJ, 953, 61. doi:10.3847/1538-4357/acdedd
  * Zhou et al. (2019) Zhou, H., Shi, X., Yuan, W., et al. 2019, Nature, 573, 83. doi:10.1038/s41586-019-1510-y
  * Zu et al. (2011) Zu, Y., Kochanek, C. S., & Peterson, B. M. 2011, ApJ, 735, 80. doi:10.1088/0004-637X/735/2/80