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

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.

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 BH (Pancoast et al. 2014;Mejía-Restrepo et al. 2018).Accurate values of M 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 spiralarm 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 sam-ples 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 distancedependent 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.

SAMPLE
We selected AGN samples with RM observations to investigate the flux ratios of Hα and Hβ emission lines, their RM results, and their joint analysis.However, simultaneous measurement of time delays for both Hα and Hβ emission lines has been realised by only a handful of AGNs observations.Fewer than 10 AGN have simultaneous velocityresolved RM measurements for both Hα and Hβ.Here we focus on the Lick AGN Monitoring Projects in 2008and 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 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α and Hβ 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 PyCALI1 , 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 ≤ 17 were chosen to ensure high-quality spectra and to determine velocityresolved 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.4meter telescope at Lijiang Observatory.These low-redshift samples allowed us to simultaneously detect Hα and Hβ 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.tract information on these broad emission lines, we have applied multi-component spectral fitting techniques, utilizing DASpec 2 , 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] λ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] λ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 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β but excluding He II λ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 2 https://github.com/PuDu-Astro/DASpecfit the host galaxy, assuming solar metallicity and an age of 11 billion years, (4) a single Gaussian for narrow lines like [O III] λλ4959,5007, and the narrow Hβ line, (5) a fourthorder Gauss-Hermite function for broad Hβ.Additionally, we fixed the widths and shifts of narrow lines to [O III] λ5007 and set the [O III] λλ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] λ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) ≥ 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γ to Hα 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α, Hβ, and Hγ (for NGC 5548, three Gaussians were separately used to fit them), but the profile of broad Hγ was bound with broad Hβ.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β, are tied to [O III] λ5007, while those of narrow emission lines with wavelengths above 6000 Å are tied to [O I] λ6300 (or [S II] λ 6718 if [O I] λ6300 is weak).Also, the flux ratios of [O III] λλ5007,4959 and [N II] λλ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] λ5007 and [O I] λ6300 (or [S II] λ 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] λ5007 widths, we convolved each spectrum to match the maximum [O III] λ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 var ), as governed by the methodology outlined in Rodríguez-Pascual et al. (1997): where S 2 is the flux variance, △ 2 is mean squared error, and ⟨F ⟩ is the mean flux over the light curve.The F var uncertainty is given by Edelson et al. (2002): where N is the total number of flux measurements within the light curve.The computed results for each light curve are compiled in Table 1.

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 τ 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 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α and Hβ 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.

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 Hα /F Hβ and F Hα /F Hγ .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α and Hβ, 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β.
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 uti- a Those five objects, whose variability characteristics prevent lag measurements from being conducted, are also illustrated in Figure 6 of Bentz et al. (2009).
lized the code JAVELIN3 , 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 ∆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α and Hβ 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α and Hβ profiles at each velocity edge, we obtained light curves of F Hα /F Hβ in these bins.Finally, we derived the weighted average velocity-resolved flux ratios by utilizing the F Hα /F Hβ light curves at each velocity edge.It is important to note that these ratios were not corrected for the time delays of Hα and Hβ 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.

Transfer Function
According to the fundamental principles of the RM technique (e.g., Blandford & McKee 1982;Welsh & Horne 1991; ) and NGC 2617.In the top panels, the black step lines represent the rms spectra of Hα.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 FHα/F Hβ values and time lags (τcent) of Hα 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α and their associated uncertainties, respectively.In the bottom panels, the relationship between FHα/F Hβ values and time lags of Hα, as obtained from various velocity bins.Bentz et al. 2010), the relationship between the light curves of the continuum and an emission line can be expressed by the following equation: where Ψ(ν, t) is the transfer function (TF) encoding the geometry and kinematics of the BLR.The terms L(ν, t) and C(t − τ ) 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 Γ-distribution (Pancoast et al. 2014): where R S = 2GM BH /c 2 is the Schwarzschild radius of a black hole with a mass M BH , R BLR is the mean BLR radius, F is the ratio of the minimum radius to the mean BLR radius, β is the shape parameter, and g is randomly drawn from a Γdistribution.For a positive variable x, the Γ-distribution is described as: where θ and α represent the scale parameter and shape parameter, respectively.In this paper, we set β =0.9, F =0.6, θ = 1, and α = 1/β 2 for simplicity.Additionally, we take into account the angular distribution of the BLR clouds, which can be described by: where θ 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 γ is considered to describe the concentration of BLR clouds (Pancoast et al. 2014).In this work, we set U γ as a random number between 0 and 1, indicating that the BLR clouds are uniformly distributed but confined within ±θ o range.
Here, we primarily model four disk-like BLRs with an opening angle of θ 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 × 10 7 M ⊙ , 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.

Anti-correlation between BD and Continuum
The BD is defined as the relative intensities of the hydrogen Balmer series emission lines.Hα and Hβ 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α and Hβ 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γ 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 ∼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 max < 0.5.Notably, our sample (including previous RM observation data) showed a significant time delay between Hα and Hβ 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β 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α and Hγ, 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 line ∝ F m con .Specifically, for Hα and Hβ, they can be expressed as powers of m1 and m2, respectively.Then the BD is determined by F Hα /F Hβ ∝ F (m1−m2) con , 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).

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 cen-LI ET AL.
-3000 0 3000 Velocity (km/s) tral 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α and Hβ 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.

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 (rmax) and the corresponding time lags for Hα and Hβ relative to the continuum.These values are denoted as rmax,Hα, r max,Hβ , τHα, and τ Hβ , respectively.
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α 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α-to-Hβ ratio at different velocities.In most cases, where gas density exceeds 10 10 cm −3 and the F Hα /F Hβ ratio falls within

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β.The results include the Spearman rank correlation coefficient (rs) and probability of the null hypothesis (p null ).
a The Spearman rank correlation analysis between continuum and BD.
b The Spearman rank correlation analysis between Hβ and BD.
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).

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 ∼ 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 ∼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β 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 BDcontinuum anti-correlation reinforces the fundamental relationship between local BD and radius, which is necessary to the velocity-resolved IM technique.

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α 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α with its higher optical depth compared to Hβ.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) ∝ 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) ∝ 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α and Hβ in BLR do exhibit differences (or different transfer functions), displaying a trend that F Hβ decreases faster than F Hα 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.

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β 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 ioniz-ing 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.

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 inflow, circularized gas, and ejecta as outflow.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 selfgravity 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 BH , as emission line profiles and time delays resulting from asymmetric BLRs significantly differ from expectations of circu-lar 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 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α and Hβ emission lines in this sample.By integrating the new technique of velocity-resolved IM with traditional velocityresolved RM, we unveiled the complex kinematics and geometry within the BLR.Our key results are as follows.
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. Tight correlations between velocity-resolved BDs and RM lags were found in 4 out of 5 AGNs with highquality 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.Furthermore, NGC 4151 exhibits an anti-correlation between velocity-resolved BDs and lags, implying complex kinematics within the BLR.
3. 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.
4. 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 radiationdriven 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.Facility: YAO:2.4mSoftware: DASpec (https://github.com/PuDu-Astro/DASpec), PyCALI (Li et al. 2014), JAVELIN (Zu et al. 2011).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 FHα/F Hβ , which is represented by the red lines.

Figure 1 .
Figure1.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α, Hβ, and flux ratios of Hα and Hβ (FHα/F Hβ ), and CCF results.Black points denote the observed data, while the orange-shaded regions illustrate the results reconstructed using JAVELIN.The FHα/F Hβ 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α, Hβ, and FHα/F Hβ , 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 FHα/F Hβ .The FHα/F Hβ 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.

Figure 3 .
Figure3.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α.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 FHα/F Hβ values and time lags (τcent) of Hα 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α and their associated uncertainties, respectively.In the bottom panels, the relationship between FHα/F Hβ values and time lags of Hα, as obtained from various velocity bins.

Figure 5 .
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 FHα/F Hβ values at different velocities and times, where FHα/F Hβ values vary with velocity and time.The black lines indicate the velocity-resolved time-averaged values of FHα/F Hβ , and the red dashed lines represent the corresponding rms values.

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

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

Figure 9 .
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.
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 ex-LI ET AL.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.

Table 1 .
Measurements of variability amplitude and relative error NOTE-The objects analyzed in this study were selected from the LAMP2008, LAMP2011, and Lijiang campaigns.The LAMP2011 datasets incorporate contributions from narrow lines.

Table 3 .
Correlation measurements 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.
LI ET AL.