Space Telescope and Optical Reverberation Mapping Project. IX. Velocity–Delay Maps for Broad Emission Lines in NGC 5548

To secure data suitable for velocity–delay mapping, NGC 5548 was the focus of an intensive monitoring campaign in 2014, the AGN Space Telescope and Optical Reverberation Mapping (AGN STORM) program. Ultraviolet (UV) spectra were obtained almost daily for 6 months with the Cosmic Origins Spectrograph on the Hubble Space Telescope (HST), securing 171 UV spectra covering rest-frame wavelengths 1130–1720 Å, including the prominent Lyα λ1216 and C iv λ1549 emission lines and the weaker Si iv λ1397 and He ii λ1640 emission lines (De Rosa et al. 2015, hereafter Paper I). During the middle two-thirds of the campaign, observations with the Swift satellite provided longer-wavelength UV, 0.3–10 keV X-ray, and optical (UBV) continuum measurements (Edelson et al. 2015, hereafter Paper II). A major ground-based campaign secured imaging photometry (Fausnaugh et al. 2016, hereafter Paper III) with sub-diurnal cadence, including the UBV and Sloan ugriz bandpasses. Optical spectroscopic observations (Pei et al. 2017, hereafter Paper V) were also obtained, with 147 spectra covering the Balmer line Hβ λ4861 and He ii λ4686.

Anomalous behavior in the emission-line response, known colloquially as the BLR Holiday, is discussed by Goad et al. (2016, hereafter Paper IV). Dehghanian et al. (2019, hereafter Paper X) present photoionization modeling using the absorption lines to diagnose how the ionizing spectral energy distribution changed during the BLR Holiday. Detailed fitting of a reverberating disk model to the HST, Swift, and optical light curves was accomplished by Starkey et al. (2016, hereafter Paper VI). The X-ray observations are discussed by Mathur et al. (2017, hereafter Paper VII). A comprehensive analysis modeling of the variable emission and absorption features is presented by Kriss et al. (2019, hereafter Paper VIII). The present manuscript, presenting velocity–delay maps derived from the spectral variations, is Paper IX.

Analysis of the STORM data sets has provided several breakthroughs and surprises that challenge our previous understanding of AGN accretion flows. One major breakthrough is the first clear measurement of interband continuum lags (Papers II and III), which can serve as a probe of the accretion disk temperature profile (Collier et al. 2001; Cackett et al. 2007). This tests a key prediction of the standard Shakura & Sunyaev (1973) disk theory, ${T}_{\mathrm{eff}}\propto {\left({M}_{\mathrm{BH}}\dot{M}\right)}^{1/4}{r}^{-3/4}$, where $\dot{M}$ is the accretion rate. The STORM results are somewhat surprising, as follows:

• 1.  
  From the continuum and broadband photometric light curves, cross-correlation analysis (Papers II and III), and detailed light-curve modeling (Paper VI), continuum lags and thus the disk size are larger than expected, by a factor of ∼3. Similarly, overlarge disks are inferred from microlensing effects in lensed quasar light curves (Poindexter et al. 2008; Morgan et al. 2010; Mosquera et al. 2013).
• 2.  
  An excess lag in the U band, which samples the Balmer continuum, suggests that the long-lag problem may be an artifact of mixing short lags from the disk with longer lags from bound-free continuum emission reverberating in the larger BLR (Lawther et al. 2018; Chelouche et al. 2019; Korista & Goad 2019). More detailed modeling is needed to see whether this hypothesis can resolve the long-lag problem and rescue the standard disk theory. A more radical proposal invokes subluminal Alfvén-speed signals that trigger local viscosity enhancements at larger radii (Sun et al. 2020).
• 3.  
  The time-delay spectrum is flatter than expected, τ ∝ λ¹ rather than τ ∝ λ^4/3. This implies a steeper temperature profile for the accretion disk, T ∝ r⁻¹ rather than r^−3/4. The best-fit power-law slope is −0.99 ± 0.03, some 7σ away from −3/4 (Paper VI). This might be evidence of nonzero stress at the innermost stable circular orbit, which can steepen the temperature profile to a slope of −7/8 (Mummery & Balbus 2020).
• 4.  
  The accretion disk spectrum, inferred from the spectrum of the variable component of the light, is much fainter than predicted using the T(r) profile inferred from τ(λ) (Paper VI). The disk surface seems to have a higher color temperature, T(r) from the time-delay spectrum τ(λ), than its brightness temperature, T(r) from the flux spectrum F_ν (λ). This low surface brightness and/or high color temperature is a further challenge to accretion disk theory. One possibility is large-grained gray dust obscuring the AGN, but that would produce a large mid-infrared excess that is not observed. Other possibilities are strong local temperature structures, or azimuthal structures in the disk thickness casting shadows on the irradiated disk surface.
• 5.  
  The light curve needed to drive continuum reverberations in the UV and optical differs in detail from the X-ray light curve (Paper VI), being smoother and lacking the rapid variations seen in the X-rays. Gardner & Done (2017) have suggested that the observations imply that the standard inner disk is largely replaced by a geometrically thick Comptonized region. Another related possibility is tilting the inner disk to align with the black hole spin (Paper VI).

These continuum reverberation results pose serious challenges to the Shakura & Sunyaev (1973) accretion disk theory, sparking new thinking on the nature of black hole accretion disks. The emission-line variations also revealed some unexpected new phenomena, as follows:

• 1.  
  There was a significant anomaly in the broad emission line behavior, the "BLR Holiday" (Paper IV). The emission lines track the continuum variations as expected in the first 1/3 of the STORM campaign, but then become fainter than expected in the latter 2/3, recovering just before the end. This anomalous period violates the expected behavior of emission lines reverberating with time delays relative to the continuum. There are also significant changes in line intensity ratios, suggesting partial covering of a structured BLR, and/or changes in the shape of the ionizing spectrum. A plausible interpretation of this BLR Holiday is that part of the BLR is temporarily obscured to our line of sight and/or shielded from the ionizing radiation by a wind outflow, launched from the inner disk, that can transition between transparent and translucent states (Dehghanian et al. 2019b).
• 2.  
  Significant broad and narrow absorption lines are seen in the UV spectra (Paper VIII). The narrow absorption lines exhibit equivalent width variations that correlate with the continuum variations. Here the time delays reflect recombination times, there being no light-travel time delays since absorption occurs only along the line of sight. The inferred density of ∼10⁵ cm⁻³ and location at ∼3 pc are compatible with clouds in the narrow-line region (NLR; Peterson et al. 2013).

The focus of this paper is an echo-mapping analysis of the emission-line variations recorded in the STORM data. Section 2 briefly describes the HST and MDM Observatory spectra and the PrepSpec analysis used to improve calibrations, and extract the mean and rms spectra and the continuum and emission-line light curves. Section 2.4 presents residuals to the PrepSpec fit, including a "Barber-Pole" pattern suggestive of a rotating structure. In Section 3, we discuss the linearized echo model and MEMEcho fit to the emission-line light curves as time-delayed echoes of the 1150 Å continuum light curve, recovering the one-dimensional delay maps Ψ(τ) for each line. To model the anomalous BLR Holiday, we extended the MEMEcho model to include slowly varying line fluxes in addition to the reverberations modeled as echoes of the driving light curve. Section 4 presents our velocity–delay maps from MEMEcho analysis of the reverberating emission-line profiles, exhibiting the clear signature of an inclined Keplerian disk with a defined outer rim and front/back asymmetry. Comparisons with previous results are discussed in Section 5, and Section 6 closes with a summary of the main conclusions.

The main results of our PrepSpec analysis are given in Figure 1 for the ultraviolet HST spectra and in Figure 2 for the optical MDM spectra, where the left column gives the mean and rms spectra and the right column gives the continuum and emission-line light curves. PrepSpec's model for spectral variations is

$$\begin{eqnarray}&&F(\lambda ,t)=A(\lambda )+B(\lambda ,t)+C(\lambda ,t),\end{eqnarray} \tag{ 1 }$$

where A(λ) is the mean spectrum, B(λ, t) models the broad emission line variations, and C(λ, t) models continuum variations. We detail these components below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

PrepSpec decomposes the mean spectrum as

$$\begin{eqnarray}&&A(\lambda )=N(\lambda )+\bar{B}(\lambda )+\bar{C}(\lambda ),\end{eqnarray} \tag{ 2 }$$

where N(λ) is the NLR spectrum, $\bar{B}(\lambda )$ is the BLR spectrum, and $\bar{C}(\lambda )$ is the continuum. These components are modeled as piecewise-cubic spline functions with different degrees of flexibility: stiff for the continuum, more flexible for the BLR, and very loose for the NLR. The emission-line components are forced to vanish outside a range of velocities around the rest wavelength of each line. After some experimentation, we set the emission-line windows to ±1500 km s⁻¹ for the NLR lines, ±10,000 km s⁻¹ for most of the BLR lines, and ±6000 km s⁻¹ for Hβ. This decomposition can be used to measure emission-line strengths, widths, and velocity profiles in the mean spectrum. However, here we use it mainly to isolate the NLR component N(λ), which PrepSpec uses to improve the flux and wavelength calibrations.

PrepSpec models the continuum variations as low-order polynomials in $\mathrm{log}(\lambda )$,

$$\begin{eqnarray}&&C(\lambda ,t)=\displaystyle \sum _{k=1}^{{N}_{c}}{C}_{k}(t)\,X{\left(\lambda \right)}^{k},\end{eqnarray} \tag{ 3 }$$

with N_c coefficients C_k (t) that depend on time. Here

$$\begin{eqnarray}&&X(\lambda )=\displaystyle \frac{\mathrm{log}({\lambda }^{2}/{\lambda }_{1}\,{\lambda }_{2})}{\mathrm{log}({\lambda }_{2}/{\lambda }_{1})}\end{eqnarray} \tag{ 4 }$$

interpolates linearly in $\mathrm{log}\lambda$ from −1 to +1 over the spectral range from λ₁ to λ₂. We adopt cubic polynomials, N_c = 4, to represent the continuum variations in the HST spectra over the rest-frame wavelengths 1130–1720 Å and linear polynomials, N_c = 2, for the optical MDM spectral range 4500–5400 Å. Lower values of N_c leave evident fit residuals, and higher values do not significantly improve the fit. PrepSpec uses the full spectral range to define continuum variations relative to the mean spectrum, rather than fitting continua to individual spectra using defined relatively line-free continuum windows.

PrepSpec models the broad emission line variations as

$$\begin{eqnarray}&&B(\lambda ,t)=\displaystyle \sum _{{\ell }=1}^{{N}_{{\ell }}}\,{B}_{{\ell }}(\lambda )\,{L}_{{\ell }}(t),\end{eqnarray} \tag{ 5 }$$

thus representing the variable component of each line ℓ as a fixed line profile B_ℓ (λ) scaled by its light curve L_ℓ (t). The light curves are normalized to $\left\langle {L}_{{\ell }}\right\rangle =0$ and $\left\langle {L}_{{\ell }}^{2}\right\rangle =1$. This constraint eliminates degeneracies between the model parameters and lets us interpret B_ℓ (λ) as the rms spectrum of the variations in line ℓ. In the same way as for the mean spectrum, the rms line profiles B_ℓ (λ) are modeled as piecewise-cubic spline functions and set to 0 outside the BLR window for that line.  This separable model for the line variations assumes that each line has a fixed line profile that varies in strength with time. Residuals to the PrepSpec fit then reveal the evidence for any changes in the line profile. Such changes contain the information we seek on the velocity–delay structure of the reverberating emission-line region and can reveal other interesting phenomena such as the rotating pattern that we discuss in Section 2.4 below.

The UV spectra are the same HST spectra discussed and analyzed with cross-correlation methods in Paper I. These spectra exhibit several narrow absorption systems that interfere with our analysis. We used the spectral modeling analysis in Paper VIII to identify wavelength regions affected by narrow absorption features and remove the narrow absorption effects. The fluxes and uncertainties in these regions are divided by the model transmission function, restoring to a good approximation the flux that would have been observed in the absence of the absorption while also expanding the error bars to appropriately reflect the lower number of detected photons.

Similarly, a Lorentzian optical depth profile provided an approximate fit to the broad wings of the geocoronal Lyα absorption. We divided the observed fluxes and their error bars by the model transmission, approximately compensating for the geocoronal Lyα absorption at moderate optical depths. The opaque core of the geocoronal line was beyond repair, and we omit those wavelengths (1214.3–1216.8 Å) from our analysis.

The main results of our PrepSpec fit to the HST spectra are shown in Figure 1, where the left and right columns show the spectral and temporal components of the model, respectively.  In Figure 1(a), the mean spectrum A(λ) is decomposed into the NLR spectrum N(λ) (orange), the BLR spectrum $\bar{B}(\lambda )$ (blue), and the continuum $\bar{C}(\lambda )$ (red). The BLR spectrum has very strong, broad Lyα and C iv emission extending to ±10,000 km s⁻¹, with weaker counterparts in Si iv and He ii. As PrepSpec failed to robustly separate Lyα and N v, we opted to model the Lyα+N v blend as a single line. The NLR spectrum is dominated by Lyα and C iv with narrow emission peaks also at N v and He ii. A few narrow absorption features remain uncorrected that will not adversely affect our analysis. Figure 1(b) shows the continuum light curves, C(λ, t), evaluated at five wavelengths across the spectrum. Continuum variations with a median absolute deviation (MAD) of 16% relative to the continuum in the mean spectrum are detected with S/N  ≈ 500. The amplitude is larger at 1130 Å on the blue end than at 1700 Å on the red end of the spectrum.  Figure 1(c) shows the rms spectra and Figure 1(d) the corresponding BLR light curves. The blue slope of the continuum variations is again evident in the rms spectrum. The BLR variations are detected with high S/N, ∼400 for Lyα, ∼300 for C iv, ∼120 for He ii, and ∼80 for Si iv. The BLR light curves generally resemble those of the continuum, but with time delays and other systematic differences that are distinct for each line. The strong Lyα and C iv lines are single peaked in the mean spectrum but double peaked in the rms spectrum, suggesting that the variations arise from a disk-like BLR. Variations are detected in N v λ1240 on the red wing of Lyα, in Si iv λ1393, and in He ii λ1640.

Optical spectra from the MDM Observatory were presented and analyzed with a cross-correlation analysis in Paper V. Ground-based spectra taken at facilities other than MDM were excluded from this analysis in order to have a consistent and homogeneous data set taken with the same instrument, same spectral resolution, and so on. The ground-based MDM spectra were taken through a 5''-wide slit and extracted with a 15'' aperture, under variable observing conditions. As a result, each spectrum has a slightly different calibration of flux, wavelength, and spectral resolution. While these residual calibration errors are most evident in the regions around narrow emission lines, they contribute to a smaller extent throughout the spectrum.

To compensate for this, the PrepSpec model M(λ, t) includes small adjustments to the calibrations:

$$\begin{eqnarray}&&M(\lambda ,t)=p(t)\left(F-{\rm{\Delta }}\lambda (t)\displaystyle \frac{\partial F}{\partial \lambda }+{\rm{\Delta }}s(t)\displaystyle \frac{{\partial }^{2}F}{\partial {\lambda }^{2}}\right).\end{eqnarray} \tag{ 6 }$$

Here the calibration-adjusted model is F(λ, t), given by Equation (1), and the small time-dependent adjustments to the calibration are parameterized by p(t) to model imperfect photometry, Δλ(t) for small changes to the wavelength scale, and Δs(t) for small changes in the spectral resolution. PrepSpec models $\mathrm{ln}p(t)$ to ensure that p(t) remains positive. The median of p(t) is set to unity, since typically a minority of the observed spectra are low owing to slit losses and imperfect pointing or variable atmospheric transparency. While PrepSpec can model $\mathrm{ln}p(t)$, Δλ(t), and Δs(t) as low-order polynomials of $\mathrm{log}\lambda$, the wavelength dependence of these calibration adjustments was not needed over the relatively short wavelength span of the MDM data analyzed here.

The main results of our PrepSpec analysis of the MDM spectra are shown in Figure 2. This optical spectral region includes the broad Hβ and He ii λ4686 emission lines and narrow [O iii] emission lines. In the mean spectrum, Figure 2(a), Hβ and [O iii] are strong but He ii is very weak. Hβ also has a narrow component. In the rms spectrum, Figure 2(b), the continuum is bluer than in the mean spectrum, the He ii line is much stronger, and the [O iii] line is very weak (if well calibrated, this emission line should not vary at all and thus should not show a signal in the rms spectrum). Note that the profile of broad Hβ is single peaked in the mean but double peaked in the rms spectrum. The continuum and emission-line light curves are shown in Figures 2(c) and (d), respectively. Maxima and minima in the He ii light curve occur a few days earlier than their counterparts in the Hβ light curve.

Figure 3 presents results of an analysis of the residuals of the PrepSpec fit to the UV HST (left) and optical MDM (right) spectra. The PrepSpec model assumes for each line a fixed line profile that changes in normalization only. The residuals to the PrepSpec fit thus present a visualization of the evidence for variations in the velocity profiles of the emission lines. They also serve as a check on the success of the absorption-line corrections, the calibration adjustments based on the narrow [O iii] emission lines, and the accuracy of the error estimates.

Zoom In Zoom Out Reset image size

The top panels, Figures 3(a) and (b), present the fitted PrepSpec model as a grayscale "trailed spectrogram," with wavelength increasing to the right and time upward. Here horizontal bands arise from the continuum variations, and vertical bands mark the locations of emission lines.  The middle panels, Figures 3(c) and (d), show residuals after subtracting the PrepSpec model from the observed spectra. There are acceptably small fine-scale residuals near the strong narrow [O iii] lines at 4959 and 5007 Å, indicating the good quality of the calibrations. In Figure 3(d) the evident patterns moving toward the center of the Hβ line arise from reverberations affecting the line wings first and then moving toward the line center. We find below that these can be interpreted as reverberation of Hβ-emitting gas with a Keplerian velocity field. There are also stationary features near 4750, 4880, and 4970 Å that decrease over the 180-day span of the observations, indicating a gradual decrease in the emission-line flux with time.

In Figure 3(c) the dominant residuals near the C iv line exhibit an intriguing helical "Barber-Pole" pattern with stripes moving from red to blue across the line profile. This Barber-Pole pattern may be present also in the Lyα residuals, but less clearly so owing to higher levels of systematic problems created by the absorption-line corrections and blending with N v. We see no clear evidence of the Barber-Pole pattern in the Hβ residuals, where the reverberation signatures are stronger. The peak-to-trough amplitude of these features in C iv is ±8% of the continuum flux density—far too large to be ascribed to calibration errors in the HST spectra.

The bottom panels, Figures 3(e) and (f), show the mean μ(λ) and rms χ(λ) of the normalized residuals, scaled by the error bars. The χ(λ) curves (blue) rise near the emission lines, where significant line profile variations are being detected, and level off in the continuum to values below unity, 0.81 for the HST and 0.89 for the MDM spectra. These low values indicate that rms residuals are smaller than expected from the nominal error bars. For the MEMEcho analysis to follow, we multiply the nominal error bar spectra by these factors.

The Barber-Pole pattern is a new phenomenon in AGNs. Manser et al. (2019) report a similar pattern of stripes moving from red to blue across the infrared Ca ii triplet line profiles arising from a thin ring or disk of gas orbiting a white dwarf. The 2 hr period is stable over several years, prompting its interpretation as due to an orbiting planetesimal perturbing the debris disk around the white dwarf.

We tentatively interpret the Barber-Pole pattern in NGC 5548 as evidence for azimuthal structure, perhaps caused by the shadows cast by the vertical structure associated with precessing spiral waves or orbiting material or streamlines near the base of a disk wind, which rotate around the black hole with a period of ∼2 yr. This 2 yr period is estimated based on the impression from Figure 3(c) that the stripes move halfway across the C iv profile during the 180-day campaign, so that 180 days is 1/4 of the period of the rotating pattern. This is clearly just a rough estimate. From the velocity–delay maps discussed in Section 4 below, we infer a black hole mass M_BH ≈ 7 × 10⁷ M_⊙ and a disk-like BLR geometry extending from 2 to 20 lt-day with an inclination i ≈ 45°. A 2 yr orbital period then occurs at R/c ≈ 4 days, or R ≈ 1000 G M_BH/c², compatible with the inner region of the BLR. The corresponding Kepler velocity is $V=\sqrt{G\,{M}_{\mathrm{BH}}/R}\approx 9000$ km s⁻¹, and this projects to $V\,\sin i\approx 7000$ km s⁻¹ for i ≈ 45°. These rough estimates are compatible with the observed velocity amplitude of the Barber-Pole stripes in the C iv residuals.

The effect must be stronger on the far side of the disk, to produce Barber-Pole features that move from red to blue across the line profile, and weaker on the near side, where they would be seen moving from blue to red. This front-to-back asymmetry might be due to a bowl-shaped BLR geometry, so that the near side of the BLR disk is strongly foreshortened. However, the velocity–delay maps discussed in Section 4 indicate that the response is stronger on the near side than on the far side of the disk. Alternatively, if the inner disk is tilted toward us, perhaps due to a misaligned black hole spin, then material orbiting there could rise above the outer-disk plane, to cast shadows on the far side of the outer disk, and then dip below the plane to avoid casting shadows on the near side of the outer disk.

Detailed modeling beyond the scope of this paper may test the viability of these and other interpretations. Further monitoring of NGC 5548 with HST may be helpful to determine whether the Barber-Pole phenomenon is stable or transient, whether its period is stable or changing, and whether the stripes always go from red to blue or sometimes from blue to red across the C iv profile.

The full spectrum of ionizing radiation is not observable, and so an observed continuum light curve, C(t), is adopted as a proxy. At each time delay τ, the responding emission-line light curve L(t) is then some nonlinear function of the continuum light curve C(t − τ) shifted to the earlier time t − τ. In addition, the observed line and continuum fluxes include constant or slowly varying background contributions from other light sources, such as narrow-line emission and starlight from the host galaxy. To model these backgrounds and account for the nonlinear BLR responses, MEMEcho employs a linearized echo model, with reference levels C₀ for the continuum and L₀ for the line flux, and a tangent-curve approximation to variations around these reference levels. Thus, the continuum light curve C(t) is decomposed as

$$\begin{eqnarray}&&C(t)={C}_{0}+{\rm{\Delta }}C(t),\end{eqnarray} \tag{ 7 }$$

and the emission-line light curve,

$$\begin{eqnarray}&&L(t)={L}_{0}+\int \,{\rm{\Psi }}(\tau )\,{\rm{\Delta }}C(t-\tau )\,d\tau ,\end{eqnarray} \tag{ 8 }$$

is a convolution of the continuum variations with a delay map Ψ(τ), giving the one-dimensional delay distribution of the emission-line response. We find that this linearized echo model fails to provide a good fit to the NGC 5548 data. We therefore generalize the model to allow a time-dependent echo background level, L₀(t). This extension is straightforward .

Maximum entropy regularization keeps the model light curves C(t) and L₀(t) and the delay maps Ψ(τ) positive and "as smooth as possible" while fitting the data. Referring to these functions generically as p(t), the entropy is

$$\begin{eqnarray}&&S(p)=\displaystyle \sum _{t}w(t)\left[p(t)-q(t)-p(t)\,\mathrm{ln}\left(p(t)/q(t)\right)\right],\end{eqnarray} \tag{ 9 }$$

measured with weights w(t) and relative to a default function q(t). We obtain q(t) by Gaussian smoothing of p(t), with an FWHM of 1, 2, and 4 days for the driving light curve, the delay map, and the echo background, respectively. These choices control the flexibility of the functions. The weights w(t) provide additional control on relative flexibility among the three functions.

For fits to reverberating spectra, the MEMEcho model simply adds a wavelength dimension to the echo light curve, L(t) → L(λ, t), to the response distribution, Ψ(τ) → Ψ(λ, τ), and to the background variations, L₀(t) → L₀(λ, t). These two-dimensional functions are then regularized with the entropy defined relative to default functions that average in both directions.

The MEMEcho fit is accomplished by iteratively adjusting the functions p to minimize

$$\begin{eqnarray}&&Q\left(p,D\right)={\chi }^{2}\left(p,D\right)-\alpha \,S\left(p\right).\end{eqnarray} \tag{ 10 }$$

Here ${\chi }^{2}\left(p,D\right)$ quantifies the "badness of fit" to the N data D, assuming Gaussian noise with known error bars. The Lagrange multiplier α controls the trade-off between fitting the data (small χ²) and being simple (large S). In practice, α is initially large, and a series of converged fits is constructed with decreasing χ² and increasing S, stopping when the fit is judged to be satisfactory (χ²/N ≈ 1) and the model not overly complex.

Figure 4 shows the results of our MEMEcho fit to five continuum and six emission-line light curves of NGC 5548. The light-curve data are from the PrepSpec analysis of the HST and MDM spectra, described in Section 2. MEMEcho fits all light curves simultaneously, recovering a model for the driving light curve C(t), and for each echo light curve a delay map Ψ(τ) and a background light curve L₀(t). The driving light curve C(t) (bottom panel of Figure 4) is the 1150 Å continuum light curve, with the reference level C₀ (red line) set at the median of the 1150 Å continuum data. Above this are 10 echo light curves (right) and corresponding delay maps (left), where the light-curve data (black points with green error bars) can be directly compared with the fitted model (blue curves). We model four continuum light curves, at 1300, 1450, and 1700 Å from the HST spectra and at 5100 Å from the MDM spectra, as echoes of the 1150 Å continuum. The reverberating emission lines are He ii λ1640 and He ii λ4686, then Hβ and Lyα, and finally Si iv and C iv. The MEMEcho fit accounts for much of the light-curve structure as echoes of the driving light curve but requires significant additional variations L₀(t) (red curves), particularly during the BLR Holiday indicated by gray shading in Figure 4.

Zoom In Zoom Out Reset image size

The fit shown in Figure 4 requires χ²/N = 1 separately for the driving light curve and for each of the echo light curves, where there are N = 171 and 147 data points for the HST and MDM light curves, respectively. The model light curves (delay maps) are evaluated on a uniform grid of times (delays) spaced by Δt = 0.5 days, linearly interpolated to the times of the observations. The delay maps span a delay range of 0–50 days.

The delay maps Ψ(τ) are of primary interest because they indicate the radial distributions from the central black hole over which the continuum and emission lines are responding to variations in the driving radiation. The continuum light curves exhibit highly correlated variations that are well fit by exponential delay distributions strongly peaked at τ = 0. The median delay, increasing with wavelength, is ∼1 day at 1700 Å and ∼5 days at 5100 Å. The echo background has only small variations, indicating that the linearized echo model is a very good approximation for the continuum light curves.

The emission-line light curves require more extended delay distributions and larger variations in their background levels. The background variations are similar, but not identical, for the six emission-line light curves. The two He ii light curves require tight delay distributions peaking at τ = 0, with half the response inside ∼5 days and 3/4 inside 10 days, and some low-level structure at 20–40 days. The background light curves have a "slow wave" with a 100-day timescale, somewhat different for the two lines, and smaller-amplitude 10-day structure. The slow-wave background for He ii λ1640 is rising from HJD 6690 to 6750 (really HJD−2,450,000), while that for He ii λ4686 is more constant. Both backgrounds then decline to minima around HJD 6800 and then rise until 6840. The constant background for He ii λ1640 prior to HJD 6690 and for He ii λ4686 after HJD 6850 is not significant since there are no data during these intervals.

The Hβ response exhibits the most extended delay distribution, with a peak at 7 days, half the response inside 14 days, 3/4 inside 23 days, minor bumps at 25 and 40 days, and falling to 0 at 50 days. The need for this extended delay map is evident in the Hβ light curve, for example, to explain the slow Hβ decline following peaks at HJD 6705 and 6745. The Lyα response is more confined than Hβ with a peak at 3 days, half the response inside 7 days, 3/4 inside 15 days, and bumps at 26 and 35 days. Si iv and C iv are similar, with peaks at 5 and 7 days, respectively.

The slow-wave backgrounds L₀(t) for all these lines fall slowly from HJD 6740 to 6820 and then rise more rapidly to a peak at HJD 6840. This corresponds approximately to the anomalous BLR Holiday period discussed in Paper IV, indicated by gray shading in Figure 4, during which the emission lines became weak relative to the continuum. Note also a smaller dip from HJD 6715 to 6740 that serves to deepen the emission-line decline between the two peaks, particularly for C iv. A small peak near HJD 6810 accounts for emission-line peaks in C iv, Si iv, and He ii λ1640 that have no clear counterpart in the continuum light curves.

Note that the model and background light curves (blue and red curves in Figure 4) exhibit numerous small spikes in addition to smoother 100-day and 10-day structure. These spikes correspond to data points that are too high or too low, relative to their error bars, to be fit by the smooth default light curve that maximizes the entropy. The largest offender is a low point in the Hβ and He ii λ4686 light curves near HJD 6837, which likely represents a calibration error. These outliers could seriously damage the delay maps. The spikes are less prominent if we relax the fit to a higher χ²/N, but then the fit to the relatively low S/N Si iv light curve is less satisfactory. Fortunately, because our model has time-dependent backgrounds that can develop sharp spikes where required, the delay maps remain relatively smooth and insensitive to these outliers.

Wavelength–delay maps Ψ(λ, τ) of the emission-line response in NGC 5548 are shown as two-dimensional false-color images in Figure 5 for the MEMEcho fit to the optical spectra from MDM and in Figure 6 for the UV spectra from HST.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the panel to the right of the 2D map, the projections Ψ(τ) give delay maps for the full wavelength range (black) and for velocity ranges centered on the rest wavelengths of the emission lines, as indicated by the colored bars above and below the map. In the panel below, the projections Ψ(λ) give the spectrum of the full response (black) and of the response in four delay ranges, 0–5 days (purple), 5–10 days (green), 10–15 days (orange), and 15–20 days (red). Velocity–delay maps centered on the six emission lines are presented in Figure 7. These two-dimensional maps show that the emission-line response inhabits the interior of a virial envelope (dashed) and exhibit structure indicating a Keplerian disk inclined by i = 45° and with an outer rim at R/c = 20 days, as discussed below. These maps and their interpretation are the main results of interest emerging from our MEMEcho analysis.

Zoom In Zoom Out Reset image size

Details from the MEMEcho fit are shown in Figures 8 and 9 for the HST and MDM spectra, respectively, with the same format as in Figure 4. The fit models as echoes the same four continuum light curves (at 1350, 1450, 1700, and 5100 Å), and now as well the continuum-subtracted emission-line variations across the full wavelength range of the HST and MDM spectra. The PrepSpec model provided the variable continuum model that was subtracted to isolate the BLR spectra used in the MEMEcho fit. As before, the proxy driving light curve C(t) is the continuum light curve at 1150 Å. The model allows echo responses over a delay range of 0–50 days and includes a time-variable background spectrum L₀(λ, t). This allows the model to account for the period of anomalous line response during the BLR Holiday (Paper IV) and any other features in the data, real or spurious, that are not easily interpretable in terms of the linearized echo model. Figure 10 presents a grayscale "trailed spectrogram" display of the slowly varying background L₀(λ, t) for the MEMEcho fit to the HST and MDM spectra. The Barber-Pole pattern of residuals and the BLR Holiday features are evident here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To regularize these fits, which include two-dimensional wavelength–delay maps Ψ(λ, τ) and varying background spectra L₀(λ, t), the entropy now steers the fit toward models with smooth spectra, as well as smooth light curves and delay maps. We actually construct a series of MEMEcho maps that fit the data at different values of χ²/N, ranging from 5 to 1. At higher χ², the fit to the data is poor and the maps are smooth. At lower χ², the fit improves and the maps develop more detailed structure. When χ² is too low, the fit becomes strained as the model strives to fit noise features (e.g., by introducing spikes in the gaps between data points in the driving light curve). The fits and maps shown here for a fit with χ²/N = 1.2 are a good compromise between noise and resolution. Our tests show that the main features interpreted here are robust to changes in the control parameters of the fit. These parameters adjust the relative "stiffness" of the driving light curve, the background light curve, the echo maps, and the aspect ratio of resolution in the velocity and delay directions.

From Figures 5 and 6, the three strongest lines—Lyα, C iv, and Hβ—have a similar velocity–delay structure, with most of their response occurring between 5 and 15 days. To first order, the response in all three lines is red–blue symmetric. This indicates that radial motions are subdominant in the BLR, as a strong inflow (outflow) component would produce shorter delays on the red (blue) side of the velocity profile (Welsh & Horne 1991).

The He ii response is largely inside 5 days and extends to ±10,000 km s⁻¹, compatible with expected radial ionization structure and virial motions. In Figure 5, we see that He ii response is broad and single peaked. He ii dominates the Hβ response in the delay slice of 0–5 days (purple), becomes subdominant at 5–10 days (green), and is almost negligible at larger delays. The He ii λ1640 and He ii λ4686 delay ranges and velocity structures are similar. There is no signature of a double-peaked structure in these two lines.

In Figure 5, the Hβ delay maps Ψ(τ) for the full profile (±6000 km s⁻¹; red) and for the line core (±1500 km s⁻¹; orange) are flat or rising from 0 to 10 days and then tail away. The Hβ response spectrum Ψ(λ) exhibits a double-peaked structure in the delay ranges of 10–15 days (orange) and 15–20 days (red), with the peaks separated by ∼5000 km s⁻¹. In the slice of 5–10 days, the Hβ response has a central peak flanked by ledges that extend to ±5000 km s⁻¹. Figure 6 shows similar double-peaked responses in Lyα and somewhat less clearly in C iv.

Clearly recognizable in the velocity–delay structure is the signature of an inclined Keplerian disk with a well-defined outer edge. The velocity–delay maps provide a plausible interpretation for the "M"-shaped variation in lag with velocity seen in the cross-correlation results (Papers I and V). The outer edges of the "M" arise from the virial envelope. The "U"-shaped interior of the "M" dips down from 20 to 5 days, and we interpret this as the lower half of an ellipse in the velocity–delay plane, which is the signature of a ring of gas orbiting the black hole at radius R = 20 lt-day.

The Hβ response exhibits the clearest signature of an ellipse in the velocity–delay plane, corresponding to an annulus in the Keplerian disk. The (stronger) near side of the annulus has a delay at $\tau =(R/c)(1-\sin i)\approx 5$ days, and the (weaker) far side extends to $\tau =(R/c)(1+\sin i)\approx 35$ days. Assuming a thin disk, the ratio gives $\sin i\approx 0.75$, or i ≈ 45°. The double-peaked velocity structure at τ ≈ 20 days then gives the black hole mass. The framework shown by dashed orange curves on Figures 8 and 9 was adjusted by eye to fit the main features. This provides rough estimates for the black hole mass M_BH ≈ 7 × 10⁷ M_⊙, the inclination i ≈ 45°, and the outer BLR radius R_out ≈ 20 lt-day. Note that the characteristic BLR response timescale, as measured by the mode or mean or median of Ψ(τ), is less than R/c at the outer edge of the BLR.

The velocity–delay structure also indicates a stronger response from the near side than from the far side of the inclined disk. We see the upper half of the ellipse only faintly in the velocity–delay map of Hβ and perhaps also for Lyα. The C iv map and the (less reliable) Si iv map show a faint response at 25–30 days that is not very clearly connected to the stronger response inside 10–15 days. If the response structure were azimuthally symmetric, the upper and lower halves of the ellipse would be more equally visible. The mean delay averaged around the ellipse would be R/c ≈ 20 days, and this is similar to typical Hβ lags seen in the past. The much shorter lags in the STORM data may be interpreted as due at least in part to an anisotropy present in 2014 that was usually much weaker or absent during previous monitoring campaigns. The near/far contrast ratio can be determined by more careful modeling.

The morphology of the velocity–delay maps indicates that the BLR in NGC 5548 is compatible with a disk-like geometry and kinematics, for a black hole mass M_BH ≈ 7 × 10⁷ M_⊙, an inclination angle i ≈ 45°, and a relatively sharp outer rim at R/c ≈ 20 days. Uncertainties in these parameters are difficult to quantify precisely because MEMEcho delivers the "smoothest positive" maps that fit the data at specified levels of χ²/N. This relatively model-independent approach does not aim to determine model-specific parameters. However, by comparing the velocity–delay maps with those of toy models, we estimate that M_BH is uncertain by ∼20%, the inclination by  ∼ 10°, and R/c by ∼10%.

Future and ongoing efforts to model the STORM data will incorporate more specific geometric and kinematic parameters and photoionization physics. These should allow for parameter estimates with better-quantified uncertainties. At this stage the value of the velocity–delay maps is to inform the dynamical modeling efforts by indicating what types of models are likely to succeed.

NGC 5548 was one of 13 AGNs monitored during the Lick AGN Monitoring Project (LAMP) 2008 RM campaign (Bentz et al. 2009; Walsh et al. 2009). Using spectra from the 3 m Shane telescope at Lick Observatory and Johnson V and B broadband photometry from a number of ground-based telescopes, the LAMP 2008 campaign secured 57 V-band continuum and 51 Hβ line epochs on NGC 5548 over a 70-day span. The centroid Hβ lag measured by cross-correlation methods is ${\tau }_{\mathrm{cent}}={4.2}_{-1.3}^{+0.9}$ days. Combined with the Hβ line width σ_line = 4270 ± 290 km s⁻¹ from the rms spectrum, the virial product is $c\,{\tau }_{\mathrm{cent}}\,{\sigma }_{\mathrm{line}}^{2}/G={14.9}_{-5.1}^{+3.7}\times {10}^{6}$ M_⊙, giving the black hole mass as ${M}_{\mathrm{BH}}={8.2}_{-2.8}^{+2.0}\times {10}^{7}\left(f/5.5\right)$ M_⊙, where f is the adopted calibration factor, uncertain by ∼0.3 dex for individual AGNs.

Pancoast et al. (2014) present results of a dynamical modeling fit to the 2008 LAMP data on NGC 5548. The caramel code employs Markov Chain Monte Carlo methods to sample the parameters of a dynamical model of the Hβ-emitting region. Taking the V-band data as a proxy for the driving light curve, caramel fits the reverberating emission-line profile variations by adjusting the spatial and kinematic distribution of Hβ-emitting clouds. This model incorporates more detailed model-specific information than the MEMEcho mapping but is flexible enough to explore a variety of inflow and outflow as well as disk-like kinematic structures. The caramel fit to NGC 5548 in 2008 finds a mean delay of just 3 days, indicating a smaller Hβ emission-line region. The caramel model fitted to the 2008 data has a thick-disk geometry, with an opening angle ${\theta }_{0}={\left({27}_{-8}^{+11}\right)}^{^\circ }$, a dominant Hβ response on the far side of the disk, and a strong inflow signature, with small delays on the red wing relative to those on the blue wing of Hβ. The inflow and far-side response in 2008 are quite different from the symmetric disk-like kinematics and near-side response evident in the MEMEcho velocity–delay maps. Despite these differences and the smaller size of the Hβ emission-line region in 2008, the inferred mass ${M}_{\mathrm{BH}}={3.9}_{-1.5}^{+2.9}\times {10}^{7}$ M_⊙ and inclination $i={\left(39\pm 12\right)}^{^\circ }$ are consistent, given their uncertainties, with the corresponding values, M_BH ≈ 7 × 10⁷ M_⊙ and i ≈ 45°, that we estimate from the structure of the Hβ velocity–delay maps from the 2014 STORM data. The mass estimate from the caramel fit may be too small owing to the Hβ response being measured relative to the optical continuum, which itself may have a delay of a few days. For example, in 2014 the V-band delay was 1.6 ± 0.5 days (Paper II), and that could boost the caramel mass estimate by ∼50%, bringing it into closer agreement with our estimate from the velocity–delay maps.

During 2015 January–July, a year after the STORM campaign, NGC 5548 was monitored with the Yunnan Faint Object Spectrograph and Camera on the 2.4 m telescope at Lijiang, China, resulting in 61 good spectra over 205 days that provide the basis for an MEM analysis yielding a velocity–delay map for Hβ (Xiao et al. 2018). This 2015 map exhibits structure remarkably similar to that seen in our 2014 map, including the virial envelope, the "M"-shaped structure with τ ≈ 10 days and V ≈ ± 2400 km s⁻¹ at the peaks of the "M," and a well-defined ellipse extending out to τ ≈ 40 days.

Xiao et al. (2018) also present Hβ velocity–delay maps constructed from 13 annual AGNWatch campaigns during 1989–2001. While several of these maps show hints of a ring-like structure, the quality of these maps is much lower owing to less intensive time sampling, making it difficult to be confident about the information that they may be able to convey.

In combination, the high-fidelity 2014 and 2015 Hβ velocity–delay maps each show a clear virial envelope and distinct ellipse centered at 20 lt-day. Thus, at both epochs the Hβ response arises from a Keplerian disk with a relatively sharp outer rim at 20 lt-day that remained stable for an interval of at least a year. In 2014 the response is weak on the top of the velocity–delay ellipse. In 2015 the response is clearly visible on the blue side and over the top of the velocity–delay ellipse and relatively weak on the red side. This indicates significant azimuthal modulation of the response that evolves, perhaps rotates, on a timescale of a year. Future velocity–delay mapping experiments to monitor this structure could be interesting to elucidate its origin and implications.

The Barber-Pole pattern uncovered here in NGC 5548 may be related to intermittent periodic phenomena seen in the subclass of AGNs that have double-peaked Balmer emission lines. In these double-peaked emitters there is a clear separation between the narrow emission lines and the broad double-peaked lines. The double-peaked velocity profiles can be modeled by emission from a Keplerian disk, with relativistic effects making the blue peak stronger and sharper than the red one and redshifting the line center. Spectroscopic monitoring of Arp 102B during 1987–1996 detected Hβ profile variations (Newman et al. 1997). In particular, during 1991–1995, the Hβ red/blue flux ratio oscillates over nearly 2 cycles of a 2.2 yr period, suggesting a patch of enhanced emission on a circular orbit within the disk. Contemporaneous monitoring of Hα during 1992–1996 shows that its red/blue flux ratio also oscillates by ±10% with a 2 yr period (Sergeev et al. 2000). A trailed spectrogram display of the Hα residuals, after subtracting scaled mean line profiles, reveals a helical Barber-Pole pattern with two stripes that move from red to blue across the double-peaked Hα profile. This is strikingly similar to what we see in the C iv profile of NGC 5548. Similar phenomena are seen in other double-peaked emitters (Gezari et al. 2007; Lewis et al. 2010; Schimoia et al. 2015, 2017). Our discovery of the Barber-Pole pattern in NGC 5548 indicates that this phenomenon is not limited to the double-peaked emitters.