A High-Quality velocity-delay map of the broad-line region in NGC 5548

NGC 5548 has been well spectroscopically monitored for reverberation mapping of the central kinematics by 19 campaigns. Using the maximum entropy method in this Letter, we build up a high-quality velocity-delay map of the H$\beta$ emission line in the light curves of the continuum and the line variations observed between 2015-2016. The map shows the response strength and lags of the velocity fields of the H$\beta$ emitting regions. The velocity-delay structure of the map is generally symmetric, with strong red and blue wings at time lag $\tau \leq 15$ days, a narrower velocity distribution at $\tau \geq 15$ days, and a deficit of response in the core. This is suggestive of a disk geometry of the broad-line region (BLR). The relatively weaker H$\beta$ response at the longer lags in the red side indicates anisotropic emission from the outer part of the BLR. We also recover the velocity-delay maps of NGC 5548 from the historical data of 13 years to investigate the long-term variability of its BLR. In general, the BLR of NGC 5548 was switching between the inflow and virialized phases in the past years. The resultant maps of seven years reveal inflow signatures and show decreasing lags, indicating that the changes in the BLR size are related to the infalling BLR gas. The other four maps show potential disk signatures which are similar to our map.


INTRODUCTION
Broad emission lines are the most prominent features of type 1 active galactic nuclei (AGNs). It is now generally accepted that these Doppler-broadened lines arise from the clouds in the broad-line region (BLR) photoionized by the central continuum radiation. Reverberation mapping (RM, see Bahcall et al. 1972;Blandford & McKee 1982;Peterson et al. 1993) has been widely used to investigate the geometry and kinematics of the BLRs in AGNs (Kaspi et al. 2000;Peterson et al. 2004;Denney et al. 2010;Bentz et al. 2010a;Grier et al. 2012;Du et al. 2018). The variations of the continuum (∆C) and the emission line (∆L) are connected as Corresponding author: Pu Du, Jin-Ming Bai dupu@ihep.ac.cn, baijinming@ynao.ac.cn the expression: where Ψ (v, τ ) is the so-called "transfer function" or velocitydelay map. The map gives the distribution of the line response over the line-of-sight velocity v and the time delay τ , and can be used to indicate the nature of the BLR. Due to the limitation of data quality, early attempts mainly focus on analyzing the time lags as a function of the line-ofsight velocity (velocity-resolved time-lag analysis, see, e.g., Barth et al. 2011a,b;Bentz et al. 2008Bentz et al. , 2009Bentz et al. , 2010aDenney et al. 2009aDenney et al. ,b, 2010Doroshenko et al. 2012;Du et al. 2016;Grier et al. 2013;Lu et al. 2016;Pei et al. 2017;De Rosa et al. 2015;Ulrich & Horne 1996). Velocity-resolved time-lag analysis has been applied to more than two dozen sources, and preliminarily reveals their BLR geometry and kinematics. However, this method measures the mean time lags in each velocity bins rather than recover the velocitydelay maps which can reveal the detailed response features. Subsequently, more advanced techniques such as the maximum entropy method (MEM, Horne et al. 1991;Horne 1994;Xiao et al. 2018), regularized linear inversion method (RLI, Krolik & Done 1995;Skielboe et al. 2015) and the dynamical modeling method (Pancoast et al. , 2014aLi et al. 2013) begin to be applied to recover the velocity-delay maps (Bentz et al. 2010b;Grier et al. 2013Grier et al. , 2017Pancoast et al. 2012Pancoast et al. , 2014bPancoast et al. , 2018Ulrich & Horne 1996;Xiao et al. 2018). Briefly speaking, MEM and RLI directly recover the velocity-delay map without adopting any specific model for the geometry and dynamics of the BLR. Dynamical modeling method fits the variations of the emission line profile by assuming specific BLR model, which can provide the velocity-delay map and the black hole mass measurement simultaneously.
As one of the most intensively observed objects in the RM study, NGC 5548 has been monitored by 19 individual campaigns, including the International AGN Watch Consortium (Peterson et al. 2002, and references therein), Bentz et al. (2007), Denney et al. (2009b), the 2008 Lick AGN Monitoring Project (LAMP2008; Bentz et al. 2009), the AGN Space Telescope and Optical Reverberation Mapping (AGN STORM;De Rosa et al. 2015;Edelson et al. 2015;Fausnaugh et al. 2016;Pei et al. 2017) and . Three velocity-resolved time-lag analyses for the Hβ line of NGC 5548 were presented, which reveal a velocity-symmetric line response and suggest that the BLR of NGC 5548 tends to be virialized or a Keplerian disk (Denney et al. 2009b;Lu et al. 2016;Pei et al. 2017). The velocityresolved time-lag measurement of LAMP2008 data shows similar time lags in each velocity bins (Bentz et al. 2009), but the corresponding dynamical modeling analysis suggests that the dynamics of the BLR is dominated by inflow (Pancoast et al. 2014b). The velocity-delay map of LAMP2008 recovered by the RLI analysis is in general consistent with the velocity-resolved result, but shows a prompt response in the red wing (Skielboe et al. 2015). In addition, Kollatschny & Zetzl (2013) investigated the BLR geometry through the profiles of the emission lines of NGC 5548, and proposed that its BLR conforms to the disk wind model (Murray & Chiang 1997;Proga & Kallman 2004), and the geometry tends to be not thick.
In 2015, we started an RM campaign of NGC 5548 to investigate its BLR physics. The mean and velocity-resolved time lags have been presented in Lu et al. (2016). This Letter recovers the velocity-delay map of NGC 5548 by using MEM. In addition, to study the long-term variations of its BLR kinematics and geometry, we also analyze the AGN watch data of 13 years and recover the corresponding velocity-delay maps.

OBSERVATIONS AND DATA REDUCTION
The spectroscopic data used in this Letter were obtained in our RM campaign during 2015 January-July, by using the Yunnan Faint Object Spectrograph and Camera mounted in the 2.4 m telescope at Lijiang Station of Yunnan Observatories of Chinese Academy of Sciences. The details of the observation and data reduction were provided in Lu et al. (2016). For completeness, we briefly introduce some general points here. (1) We oriented the long slit (2 .5 wide) to take the spectra of NGC 5548 and a nearby comparison star simultaneously, and calibrate the spectra of the object by using the comparison star. This procedure gives a calibration accuracy of ∼ 2% for NGC 5548 (see more details in ); (2) The instrument broadening is ∼ 500 km s −1 (see Du et al. 2016), which is much less than the FWHM of the Hβ line in the mean spectrum (∼ 9912±362 km s −1 , see Table 4 in Lu et al. 2016); (3) The time span of this observation is 205 days (with 61 epochs); (4) We used a spectral fitting scheme to measure the continuum and Hβ line fluxes (see details in Hu et al. 2015;Lu et al. 2016), which allows us to remove the irrelevant spectral components as the He II and Fe II emission, the narrow lines, and the host galaxy. In the MEM analysis of our observation, we use the broad Hβ profiles after subtracting those irrelevant components.

MEM FITTING
MEM is an effective method for fitting data without assuming any algebraic form of the model. The method has been applied to the RM observations, and successfully recovered their corresponding velocity-delay maps of about 11 sources (see Ulrich & Horne 1996;Bentz et al. 2010b;Grier et al. 2013;Xiao et al. 2018). The principle and equations of the MEM have been described in Horne (1994) and Xiao et al. (2018), and we make a brief overview here. Generally speaking, MEM introduces a linearized echo model (Equation 1; for convenience, we discretize it) (2) to fit the observed continuum light curve and the variations of the emission-line profiles. The velocity-delay map Ψ (v i , τ j ), the continuum C(t k ), and the background spectrumL 0 (v i ) are treated as parameterized model in the MEM fitting. Herē L 0 (v i ) is designed to account for the non-variable component of the emission line, andC 0 is a reference continuum level, where we adopt the median of the continuum flux asC 0 . The MEM fitting is accomplished by minimizing where  the data D m (D m and M m ( p) include the emission line and continuum light curves), S = n [p n − q n − p n ln(p n /q n )] is the entropy which controls the "simplicity-of-modeling" by minimizing the differences between the model parameters p n and the "default image" q n . Here m and n denote the numbers of the observational points and the parameters in the model, respectively. The parameter α controls a trade-off between the "goodness-of-fitting" and the "simplicity-ofmodeling", which means increasing α smooths the MEM model and leads to larger χ 2 /N , and vice versa. The model parameter p n includes Ψ (v i , τ k ), C(t k ) andL 0 (v i ), and q n is designed as the geometric mean of p n . For one-dimensional model components (C(t) andL 0 (v)), we define where x is t or v for C(t) orL 0 (v), respectively, and for two-dimensional model (Ψ (v, τ )): Here A is a parameter which assigns the weight and controls the aspect ratio of Ψ (v i , τ k ) in v and τ direction. Increasing A smears out the fine structures along the v direction, and vice versa. In this way, the total entropy can be written as where W Ψ and W C are weight parameters which control the relative "stiffness" of L(v i , t k ), C(t k ) andL 0 (v i ). In the MEM fitting, α, A, W Ψ and W C are the user-controlled parameters. The selection of these parameters has been discussed in Xiao et al. (2018). Similar to Grier et al. (2013) and Xiao et al. (2018), we first model the continuum by using the damped random walk (DRW, Li et al. 2013;Zu et al. 2013), then use the resulting highly sampled continuum in the MEM fitting instead of the original one. One benefit of this approach is that the DRW model can give reliable uncertainties. Additionally, the DRW model can be used to extrapolate the continuum light curve to times shortly before the campaign began, and thus provide better constraints on the MEM continuum modeling.
In Figure 1, we demonstrate the MEM fitting of the emission-line light curves at some uniformly-spaced (in velocity space) wavelengths, and draw the corresponding one-dimensional transfer functions in the left panels. The reconstruction of the continuum is shown in the bottom panel. It is obvious that the transfer functions of the line wings exhibit simple structures (basically show only one dominant peak), whereas the transfer functions around the line core are relatively complex with at least two peaks.
In order to better illustrate the overall fitting, we compare the time series of the line profiles and the corresponding MEM recovery in Figure 2. In general, the model fits nicely with the Hβ line profiles at all epochs, and the χ 2 /N of the overall fitting is 1.296. In addition, the residual (bottom panel of Figure 2) shows some weak signals (at ∼ 4959Å and ∼ 4861Å), which are coming from the imperfect [O III] λ4959 and Hβ narrow line subtractions. Figure 3 plots the Hβ velocity-delay map of NGC 5548 recovered from our data. We mark the wavelengths, which are selected to be shown in Figure 1, in its top axis. The map shows a symmetric "bell" shape, with a wide velocity distribution at short lags ( 15 days) and a narrower velocity dispersion at longer lags. The broad wings extend to ∼ ±7200 km s −1 , while the response in the line core extends to ∼ 48 days. In addition, there is a hollow in the core of the map, and the response is relatively weaker at ∼ [2400 km s −1 , 28 days]. For comparison, we plot the "virial envelope" v 2 = GM • /cτ in Figure 3 with dotted lines, where M • is the black hole mass, G is the gravitational constant, and c is the speed of light. Here we adopt M • = 8.71 × 10 7 M in Lu et al. (2016). The line re-  sponse is compatible with the envelope, despite the weak "spine" at [v ∼ 6000 km s −1 , τ ∼ 12 − 30 days] and the "blob" close to the line core at ∼ 60 days which are affected by the residuals coming from [O III] λ4959 and Hβ narrow line subtractions, respectively. There is a weak response at ∼ [−3000 km s −1 , 50 days], by comparing with the error map in Figure 5 (see the next section), this weak feature is significant. Its origin and evolution merit further investigations.

THE VELOCITY-DELAY MAP
In Figure 3, we also demonstrate a comparison between the mean time lags at different velocities in our map and the velocity-resolved time-lag measurements. Here the velocityresolved time lags are derived by dividing the Hβ profile into 33 uniformly-spaced (450 km s −1 ) bins, and crosscorrelating the Hβ light curve in each bin with the continuum (see more details in Lu et al. 2016). This result is essentially identical to that presented in Lu et al. (2016), although the Hβ profile is divided into narrower bins. The corresponding maximum cross-correlation coefficients (r max ) are marked by the gray-scale bars. To do the comparison, we convolve the velocity-delay map with the autocorrelation function (ACF) of the continuum (the output is identical to the CCF), and calculate the centroid around the peaks  (> 80%) of the outputs as the mean time lags. As expected, in general, the two results are consistent. The double-peaked structure of the velocity-resolved lags is similar to what was found by Pei et al. (2017) for the 2014 campaign.
Theoretically, a virialized BLR produces a velocitysymmetric signature like "bell" shape. This is because the velocity-delay structure of a Keplerian orbit is an ellipse, the orbits at inner (outer) radii of a virialized BLR produce ellipse structures on the map with wider (narrower) velocity distribution at shorter (longer) time delay, and the map is confined within the "virial envelope" (e.g., Bentz et al. 2010b;Grier et al. 2013;Xiao et al. 2018). In particular, the map of an inclined Keplerian disk has a lack of response in the core, and is different from the signature of a spherical shell, which has a filled 'bell' shape (e.g., see Figure 1 of Horne et al. 2004 and Figure 14 of Grier et al. 2013).
The "bell-like" velocity-delay map of NGC 5548 implies that its BLR is probably an inclined disk. It is unlikely to explain the map as a spherical shell geometry, because of the response deficit from ∼ 20 days to ∼ 32 days. The response of the map on the red side at ∼ [2400 km s −1 , 28 days] is relatively weaker. Considering that the response of this area comes from the outer part of the BLR, such evidence indicates that the Hβ response at the outer radius may be more anisotropic or inhomogeneous.

LONG-TERM VARIABILITY
In order to investigate the long-term variability of the BLR kinematics and geometry, we compile the historical spectroscopic data from the AGN Watch archive 1 , which is by far the largest optical monitoring project of NGC 5548. This project involves 13 observing campaigns from Dec 1988 to Sept 2001, and each campaign has a time span of more than ∼ 280 days. The spectra are calibrated by using the [O III] narrow emission line as in Peterson et al. (2002). We did not apply the fitting scheme described in Section 2 to the AGN Watch data.
In Figure 4, we present the 13 velocity-delay maps recovered from the AGN Watch data together with our map (the same as in Figure 3). For comparison, we draw a simulated velocity-delay map for an inclined Keplerian disk in the bottom right panel of Figure 4. The disk is inclined (i = 45 • to the observer) with an inner radius of R in = 3 lt-days and an outer radius of R out = 28 lt-days, the emissivity distribution of the BLR clouds is assumed to be ∝ R −1 . The AGN Watch maps are denoted as AGN Watch 01-13. We apply the flux randomization (see details in Peterson et al. 1998), which modify the flux of each datum by a random Gaussian deviate within the flux uncertainty, and use the Monte Carlo (MC) simulations to calculate the uncertainties of the 1 http://www.astronomy.ohio-state.edu/ ∼ agnwatch/ velocity-delay maps. The error maps are shown in Figure 5. Given the error maps, the differences between the velocitydelay maps are significant. From the AGN Watch 01 to 05, the maps show continuing inflow signatures with longer lags at the blue end and shorter lags toward the red end. Interestingly, at this period, the average time lags are generally decreasing as well (see Table 8 in Peterson et al. 2002). The AGN Watch 07 to 08 are also dominated by inflow signatures with a decrease in the time lags. It implies that the shrink of its BLR may correlate to the inflow dynamics. It has been illustrated that the BLR size of NGC 5548 follows its continuum luminosity (Kilerci Eser et al. 2015). However, subsequent study reveals that the variation of the BLR size lags 2.35 +3.47 −1.25 yrs behind the luminosity change, and this lag is similar to the dynamical timescale (∼ 2.1 yrs) of its BLR . The recombination timescale is insufficient to explain the variation of the BLR time lag, the change of its BLR kinematics is needed . The maps of the AGN Watch 06, 09, 11, and 13 generally show symmetric signatures with a paucity of response in the cores, which are similar to the signature of an inclined disk found in our map. They reveal that the BLR is disk-dominated in these periods. However, the detections of the response are limited by the data quality. The rest of the maps are not well recovered due to the low sampling rate or the small Hβ variability (see Table 6 and 7 in Peterson et al. 2002). Figure 4 shows a chronological series of the velocity-delay maps, indicating transitions between inflows and virialized status. This provides evidence for the BLR origin from the tidally disrupted clumps from the torus (Wang et al. 2017).

SUMMARY
We present the high-quality Hβ velocity-delay map of NGC 5548 recovered from our RM campaign in 2015. The map clearly shows a symmetric "bell-like" signature with a lack of response in the core. Such a structure is in accord with the predicted map of a Keplerian disk. The weaker response in the red than in the blue side at ∼ 28 days of the map indicates that the response at the outer radius of the BLR may be anisotropic or inhomogeneous. We also show the velocitydelay maps constructed from the 13-years AGN Watch data. The maps of the seven years reveal that the decreasing BLR size is probably related to the inflowing BLR gas. The other four maps show potential disk signatures which are consistent with our map. The velocity-delay maps of NGC 5548 imply that its BLR was switching between the inflow and virialized status in the past years.
We thank the anonymous referee for constructive suggestions. We thank Hong-Tao Liu for the supports of providing the computing resources, and Michael S. Brotherton for his helpful suggestions that improved the manuscript. We acknowledge the support of the staff of the Lijiang 2.4m telescope. Funding for the telescope has been provided by CAS and the People's Government of Yunnan Province. This research is supported by the National Key R&D Program of China (grants 2016YFA0400701 and 2016YFA0400702), by NSFC through grants NSFC-11503026, -11233003, -11573026, -11703077, -11773029, by