Space Telescope and Optical Reverberation Mapping Project.III. Optical Continuum Emission and Broad-Band Time Delays in NGC 5548

Space Telescope and Optical Reverberation Mapping Project.
III. Optical Continuum Emission and Broad-Band Time Delays in NGC 5548

M. M. Fausnaugh1    K. D. Denney1,2,3    A. J. Barth4    M. C. Bentz5    M. C. Bottorff6    M. T. Carini7    K. V. Croxall1,2    G. De Rosa1,2,8    M. R. Goad9    Keith Horne10    M. D. Joner11    S. Kaspi12,13    M. Kim14    S. A. Klimanov15    C. S. Kochanek1,2    D. C. Leonard16    H. Netzer12    B. M. Peterson1,2    K. Schnülle17    S. G. Sergeev18    M. Vestergaard19,20    W.-K. Zheng21    Y. Zu1,22    M. D. Anderson5    P. Arévalo23    C. Bazhaw5    G. A. Borman18    T. A. Boroson24    W. N. Brandt25,26,27    A. A. Breeveld28    B. J. Brewer29    E. M. Cackett30    D. M. Crenshaw5    E. Dalla Bontà31,32    A. De Lorenzo-Cáceres10    M. Dietrich33,34    R. Edelson35    N. V. Efimova15    J. Ely8    P. A. Evans9    A. V. Filippenko21    K. Flatland16    N. Gehrels36    S. Geier37,38,39    J. M. Gelbord40,41    L. Gonzalez16    V. Gorjian42    C. J. Grier    1,25,26    D. Grupe43    P. B. Hall44    S. Hicks7    D. Horenstein5    T. Hutchison6    M. Im45    J. J. Jensen19    J. Jones5    J. Kaastra46,47,48    B. C. Kelly49    J. A. Kennea24    S. C. Kim14    K. T. Korista50    G. A. Kriss8,51    J. C. Lee14    P. Lira52    F. MacInnis6    E. R. Manne-Nicholas5    S. Mathur1,2    I. M. MHardy53    C. Montouri54    R. Musso6    S. V. Nazarov18    R. P. Norris5    J. A. Nousek25    D. N. Okhmat18    A. Pancoast55,56    I. Papadakis57,58    J. R. Parks5    L. Pei4    R. W. Pogge1,2    J.-U. Pott17    S. E. Rafter13,59    H.-W. Rix17    D. A. Saylor5    J. S. Schimoia1,60    M. Siegel24    M. Spencer11    D. Starkey10    H.-I. Sung14    K. G. Teems5    T. Treu49,61,62    C. S. Turner5    P. Uttley63    C. Villforth64    Y. Weiss13 J.-H. Woo45    H. Yan65    and S. Young35

We present ground-based optical photometric monitoring data for NGC 5548, part of an extended multi-wavelength reverberation mapping campaign. The light curves have nearly daily cadence from 2014 January to July in nine filters (BVRI and ugriz). Combined with ultraviolet data from the Hubble Space Telescope and Swift, we confirm significant time delays between the continuum bands as a function of wavelength, extending the wavelength coverage from 1158 Å to the band ( Å). We find that the lags at wavelengths longer than the V band are equal to or greater than the lags of high-ionization-state emission lines (such as He ii and ), suggesting that the continuum-emitting source is of a physical size comparable to the inner broad-line region (BLR). The trend of lag with wavelength is broadly consistent with the prediction for continuum reprocessing by an accretion disk with . However, the lags also imply a disk radius that is 3 times larger than the prediction from standard thin-disk theory, assuming that the bolometric luminosity is 10% of the Eddington luminosity (). Using optical spectra from the Large Binocular Telescope, we estimate the bias of the interband continuum lags due to BLR emission observed in the filters. We find that the bias for filters with high levels of BLR contamination () can be important for the shortest continuum lags, and likely has a significant impact on the u and U bands owing to Balmer continuum emission.

galaxies: active — galaxies: individual (NGC 5548) — galaxies: nuclei — galaxies: Seyfert

1 Introduction

The continuum emission of radio-quiet active galactic nuclei (AGN) is believed to originate in an accretion disk around a supermassive black hole (SMBH). At accretion rates and masses appropriate for SMBHs, geometrically thin, optically thick accretion disks have maximum temperatures of K, naturally accounting for the characteristic peak of AGN spectral energy distributions (SEDs) in the far ultraviolet (UV; Burbidge, 1967; Shakura & Sunyaev, 1973; Shields, 1978). However, a large variety of competing models of the accretion flow exist, such as thick-disk geometries at extremely super- or sub-Eddington accretion rates (Abramowicz et al., 1988; Narayan & Yi, 1995). In addition, AGN exhibit nonthermal X-ray emission, which requires a hot plasma component or “corona” (e.g., Haardt & Maraschi, 1991; Chakrabarti & Titarchuk, 1995). The potential configurations and complex interplay between the hot corona and accretion disk admit a wide range of models with many free parameters, and searching for the unique observational signatures of a given disk model is very challenging (Sun & Malkan, 1989; Laor et al., 1997; Koratkar & Blaes, 1999; Vestergaard & Wilkes, 2001; Telfer et al., 2002; Kishimoto et al., 2004, and references therein).

Reverberation mapping (RM; Blandford & McKee, 1982; Peterson, 1993, 2014) can provide insight into the structure of the accretion disk, and has become a standard tool for AGN astrophysics over the last 25 years (Clavel et al., 1991; Peterson et al., 1991; Horne et al., 1991; Kaspi et al., 2000; Peterson et al., 2004; Bentz et al., 2009; Denney et al., 2010; Grier et al., 2013b; Pancoast et al., 2014; Pei et al., 2014; Barth et al., 2015, and references therein). The basic principle of RM is that emission at two different wavelengths is causally connected, so that the time delay (or lag) between two light curves represents the light-crossing time within the system, and thereby provides a straightforward measurement of the system’s physical size. For example, because the AGN continuum powers the prominent emission lines observed in Seyfert galaxy/quasar spectra, the time delays between continuum and broad-line light curves are commonly used to determine the physical extent of the line-emitting gas (the so-called broad-line region, BLR).

In a similar way, RM techniques can be used to constrain the physical processes governing AGN continuum emission. X-ray emission from the corona may irradiate and heat the accretion disk. If the corona is relatively compact and centrally located, the UV and optical emission would be expected to respond to the incident X-ray flux, “echoing” the X-ray light curve after a time delay corresponding to the light-travel time across the disk (Krolik et al., 1991). On the other hand, X-ray light curves would be expected to lag behind UV and optical light curves if the X-rays are produced by Comptonization of thermal UV/optical disk photons (Haardt & Maraschi, 1991). Observational investigations of the relation between X-ray and UV/optical emission have produced ambiguous results. X-rays have been found to lead the optical emission by one to several days in some objects (e.g., Arévalo et al., 2009; Breedt et al., 2010; Shappee et al., 2014; Troyer et al., 2016), but the X-ray variability on long ( year) timescales cannot always account for the optical variations (Uttley et al., 2003; Breedt et al., 2009). In addition, other studies find no long-term X-ray/optical correlations (Maoz et al., 2002), or find optical variations that lead the X-rays on shorter timescales ( days, Marshall et al., 2008).

RM can also reveal information about the size and structure of the continuum-emitting source. Emission from different portions of the disk peaks at different wavelengths depending on the local disk temperature. By translating the wavelength of the continuum emission into a characteristic temperature, time delays between continuum light curves can be used to map the temperature profile of the disk. The first statistically significant interband time delays were found in NGC 7469 by Wanders et al. (1997) and Collier et al. (1998). Sergeev et al. (2005) carried out intensive optical monitoring of 14 AGN and found evidence that longer wavelengths lag shorter-wavelength emission. More recent continuum RM campaigns have used the Swift observatory (Gehrels et al., 2004) to obtain unprecedentedly well-sampled light curves across X-ray, near-UV, and optical wavelengths: Shappee et al. (2014) observed NGC 2617 with Swift on a nearly daily basis for several months in 2014, while McHardy et al. (2014) monitored NGC 5548 with day cadence for approximately 2 years (excepting seasonal gaps). These studies found trends of lag with wavelength that are well fit by the expectation for X-ray/far-UV reprocessing.

The present study is the third in a series describing the results of the AGN Space Telescope and Optical Reverberation Mapping (STORM) project, an intensive, multi-wavelength monitoring campaign of NGC 5548. The AGN STORM campaign is anchored by daily far-UV observations using the Cosmic Origins Spectrograph (COS; Green et al. 2012) on the Hubble Space Telescope (HST). De Rosa et al. (2015, hereinafter Paper I) give a complete introduction to the project and an analysis of the HST data. The COS program was complemented by a four-month broad-band photometric monitoring campaign using Swift, the first results of which are presented by Edelson et al. (2015, hereinafter Paper II). The Swift campaign achieved day cadence and detected significant lags between the UV and optical continua, which follow the expected lag-wavelength relation of a thin accretion disk ().

Supplementing these space-based observations are ground-based optical monitoring programs. The present study details the optical broad-band photometric monitoring component, which extends the analysis in Paper II using data in nine optical filters with day cadence for seven months. The similarly intensive ground-based spectroscopic monitoring will be presented by Pei et al. (in prep., hereinafter Paper V). In terms of cadence, temporal baseline, and wavelength coverage, the combination of UV and optical observations of the AGN STORM project represents the most complete RM experiment ever conducted.

The present work has three primary goals. The first is to directly compare the far-UV and optical light curves of NGC 5548 over a concurrent monitoring period. The far-UV light curve (1350 Å) is expected to closely trace the true ionizing continuum (912 Å), while the optical continuum (5100 Å) appears to be delayed and somewhat smoothed compared to the UV emission. Since ground-based RM campaigns use the optical continuum as a proxy for the driving continuum light curve, understanding how the continuum emission changes as a function of wavelength is important for understanding any systematic effects in optical RM experiments. The second goal is to search for time delays between the UV and optical data, in an attempt to probe the structure of the continuum-emitting region. However, because broad-band filters pick up spectral features that arise in the BLR (e.g., strong emission lines), and these features have large lags relative to the underlying continuum (several days for a Seyfert galaxy such as NGC 5548), interband lags estimated from broad-band photometry may be biased indicators of the accretion-disk size. Therefore, our final goal is to estimate the impact of BLR emission on the observed interband time delays.

In §2, we describe the observations, data reduction, flux calibration, and general properties of the ground-based photometric light curves. In §3, we describe our time-series analysis, measuring the lag as a function of wavelength of the broad-band filters. We then explore the impact of BLR emission on the observed interband lags in §4. Finally, in §5, we discuss our results, and we summarize our conclusions in §6. Where relevant, we assume a standard cosmological model with , , and km s (Komatsu et al., 2011).

2 Observations and Data Reduction

In conjunction with the HST COS UV RM campaign described in Paper I, NGC 5548 was observed between 2013 December and 2014 August by 16 ground-based observatories in optical broad-band filters: Johnson/Cousins BVRI and Sloan Digital Sky Survey (SDSS) ugriz. A short description of each telescope, the relevant imager, and the number of contributed epochs is given in Table LABEL:tab:telescopes. All observatories followed a common reduction protocol: images were first overscan-corrected, bias-subtracted, and flat-fielded following standard procedures. The reduced data, as well as nightly calibration frames, raw images, and observing logs, were then uploaded to a central repository, and the image quality was assessed by eye. Images taken in reasonable atmospheric conditions and free of obvious reduction errors were analyzed as described below.

Table 1: Bohyunsan Optical Astronomy BOAO 1.8m e2v CCD231-84 021 March-April V 5 aaaObservatory Crimean Astrophysical CrAO 0.7m AP7p CCD 176 Dec-June BVRI 76 aaaObservatory Fountainwood Observatory FWO 0.4m SBIG 8300M 035 Jan-August V 60 Hard Labor Creek Observatory HLCO 0.5m Apogee USB/Net 075 April-June V 27 La Silla Observatory GROND 2.2m Gamma-ray Burst Optical/ 033 Feb-July griz 6 Near-IR Detector Las Cumbres Observatory LCOGT 1.0m SBIGSTX-16803/ 023 Jan-August BV ugriz 263 aaaGlobal Telescope Network Sinistro CCD-486BI 039 Lick Observatory Katzman KAIT 0.8m AP7 CCD 080 Jan-July V 80 aaaAutomatic Imaging Telescope Liverpool Telescope LT 2.0m e2v CCD 231 015 Feb-August ugriz 120 Maidanak Observatory MO15 1.5m SNUCAM 027 April-August BVR 45 Mt. Laguna Observatory MLO 1.0m CCD2005 041 June-August V 10 Mt. Lemmon Optical LOAO 1.0m KAF-4301E 068 Feb-July V 26 aaaAstronomy Observatory Nordic Optical Telescope NOT 2.5m e2V CCD42-40 019 April V 3 Robotically Controlled RCT 1.3m SITe CCD 030 Dec-May BV 55 aaaTelescope Svetloe Observatory SvO 0.4m ST-7XME CCD 200 Jan-May BVRI 49 West Mountain Observatory WMO 0.9m Finger Lakes PL-3041-UV 061 Jan-July BVR 44 Wise Observatory WC18 0.5m STL-6303E CCD 147 Dec-July BVRI 126 \@ifx@empty\@ifx@empty

2.1 Differential Photometry

The analysis is based on the ISIS image-subtraction software package (Alard & Lupton, 1998). Images are first registered to a common coordinate system, and the images with the lowest backgrounds and best seeing are combined into a high signal-to-noise ratio (S/N) “reference” image. The other images are then rescaled to match the effective exposure time of the reference image. Next, the reference image is convolved with a spatially variable kernel to match the point-spread function (PSF) of each individual epoch, and then subtracted to leave the variable flux in each pixel. We use ISIS’s built-in photometry package to extract light curves from the subtracted images at the position of the AGN in NGC 5548, in units of differential counts relative to the reference image. Because each telescope/filter/detector combination has slightly different properties (pixel scales, fields of view, gains, etc.), we built reference frames and subtracted images for each unique dataset. This procedure corrects for variable seeing conditions and removes nonvariable sources such as host-galaxy starlight, allowing a clean measurement of the variable AGN flux.

2.2 Measurement Uncertainties

The formal errors found by ISIS sometimes underestimate the full uncertainties because they only account for local Poisson error contributions. In order to estimate more reliable measurement uncertainties, we examined the residual fluxes of stars in the subtracted images, and rescaled the formal ISIS errors to be consistent with the observed scatter. Our method is similar to that of Hartman et al. (2004, §4.1).

We first used ISIS to extract differential light curves at the positions of each unsaturated star in the reference images. For stars with constant flux in time, the distribution of residual fluxes at each epoch serves as an estimate of the uncertainty in the subtraction. Since we are only concerned with the magnitude of the residuals, we first take their absolute value. We then divide these values by their formal ISIS uncertainties, so that the resulting ratios indicate the factor by which the true uncertainties are underestimated. We set a minimum value of 1.0 for this ratio, since the uncertainty cannot be smaller than the local photon noise. Finally, we multiply the formal uncertainty for the AGN at the matching epoch by the median of the rescaling factors from all stars. The procedure ensures that the measurement uncertainty in a given image is consistent with the observed scatter of the subtracted stars. The median rescaling factor for all images was 2.9, while 75% of the rescaling factors are less than 6.6 and 98% are less than 25.0. The remaining 2% have rescaling factors between 30 and 87. The poorest subtractions result when ISIS cannot accurately construct the image PSF, usually because the image has too few stars.

To assess the effectiveness of this method, we adjusted the stars’ uncertainties by the derived rescaling factor for each image, and then checked the goodness-of-fit for a constant-flux model of each star. The goodness-of-fit is calculated by


where is the number of degrees of freedom of the fit, is the counts in the light curve at epoch , is the rescaled uncertainty, and is the mean counts of the light curve. 90% of the rescaled values of are between 0.32 and 2.09, and the distribution peaks at , somewhat lower than would be expected for purely Gaussian statistics. This may indicate that our rescaling method is slightly overestimating the measurement uncertainties. However, given our large dataset, we can afford to be conservative in this regard.

Data from the Katzman Automatic Imaging Telescope (KAIT; Filippenko et al., 2001) and the u-band data from the Liverpool Telescope (LT; Steele et al., 2004) required a different treatment since these images have 10 stars or fewer, which are not enough to provide robust estimates of the error-rescaling factors. We instead calculated global rescaling factors from all available epochs, rather than individual corrections from single images. Using Equation 1, we calculate for all available stars, using the unscaled ISIS uncertainties for . We then multiplied the uncertainties of the AGN light curve by the average value of . We found that the mean rescaling factor was 8.99 for the KAIT data and 2.23 for the u-band LT data. Although this method does not account for epochs with high-quality subtractions, we find the cautious approach preferable to underestimating the uncertainties.

2.3 Intercalibration of Light Curves

In order to combine differential light curves in the same filter but from different telescopes, it is necessary to intercalibrate the light curves to a common flux scale. This accounts for the different mean flux levels and analog-to-digital unit (ADU) definitions between the reference images, as well as small differences in filter transmission functions, detector efficiencies, etc., of the many telescopes. We model the difference of any two light curves by a multiplicative rescaling factor and an additive shift . While it is trivial to solve for these parameters by matching epochs where the fluxes are known to be equal, no two observations occur at precisely the same time and it is therefore necessary to interpolate the light curves. Furthermore, this method can only treat two light curves at a time, and therefore loses information by ignoring the global probability of the ensemble calibration parameters for all telescopes. In order to address both of these problems, we built a full statistical model of the intercalibrated light curve using the software package JAVELIN, following the SPEAR formalism of Zu et al. (2011).

JAVELIN models the light curves as a damped random walk (DRW). Although recent studies using Kepler light curves have shown that the DRW overpredicts the amplitude of AGN continuum variability on short timescales (Edelson et al., 2014; Kasliwal et al., 2015), the DRW provides an adequate description of the observed light curves for the noise properties and cadence/timescales of this study (we quantitatively verify this claim in the Appendix, but see also Kelly et al. 2009; MacLeod et al. 2010; Zu et al. 2013). In brief, points sampled from a DRW have an exponential covariance matrix, which is described by an amplitude that characterizes the strength of short-term variations, and a damping timescale over which the light curve becomes decoherent. Using a Markov Chain Monte Carlo (MCMC) calculation, we simultaneously fit for the shifts and rescaling factors of all light curves in a single filter. We also fit for , but our light curves do not have a sufficiently long temporal baseline to meaningfully constrain . We therefore fixed = 164 days, so as to match the value determined from multiyear historical light curves of NGC 5548 (Zu et al., 2011). The model provides a well-defined and self-consistent means of interpolating all the light curves simultaneously (see Zu et al. 2011; Li et al. 2013 for further details).

Our fitting procedure requires one light curve to be chosen to define the flux scale and mean flux level of the resulting intercalibration, so that this reference light curve is assigned a shift of and a rescaling factor of . In the Johnson BVRI bands, we use the Wise C18 (WC18; Brosch et al., 2008) data as the calibration light curve, owing to its dense temporal sampling, long baseline, and large number of comparison stars (400). For the SDSS ugriz bands, we use the LT light curves, since they have the longest baseline and most complete time sampling.

Uncertainty in the intercalibration parameters for a given telescope contributes to the final measurement uncertainty. For a flux measurement at epoch , the calibrated measurement is , and standard error propagation shows that the uncertainty introduced per point is . Since and are usually anticorrelated, is often small compared to the uncertainties from image subtraction. However, this is not always the case for telescopes with very small numbers of observations, so we calculated , , and from the posterior distributions of these parameters for each telescope, and added in quadrature to the rescaled ISIS uncertainties for each epoch. This treatment is very conservative, since the intercalibration uncertainty is strongly correlated between points from a single telescope. A summary of the mean intercalibration uncertainties is given in Table LABEL:tab:caldat.

The choice of reference light curves defines the physical flux level of the AGN from the corresponding ISIS reference image (WC18 and LT). We convert the intercalibrated differential light curves to physical flux units by performing aperture photometry on these reference images. For the AGN and all unsaturated stars in the field, we measured the flux enclosed in a 50 radius circular aperture, and converted the summed fluxes to instrumental magnitudes. This means that the host-galaxy light within the aperture is included in the measurement of the AGN flux, and this issue is discussed in §2.4. The background sky level was estimated from an annulus with inner/outer radius of 14/29 for the stars, and 118/132 for the AGN (so as to avoid light from the host galaxy).

We then matched all stars to the SDSS Data Release 7 catalog (Abazajian et al., 2009), and computed the offset between instrumental magnitudes and the SDSS AB magnitudes. We did not find any significant color terms in the flux calibration from the comparison stars, although the spectral slope of the AGN would be a poor match to such color terms regardless of their small values. For the Johnson/Cousins BVRI bands, we determined the comparison-star magnitudes using the filter-system transformations given by Fukugita et al. (1996), and converted these to AB magnitudes using Fukugita et al. (1996) Table 8. The filter transforms have an uncertainty of mag, which we adopt as a floor for the BVRI flux-calibration uncertainty. The final flux-calibrated light curves are shown in Figure 1 and given in Table LABEL:tab:lightcurves.

Figure 1. : BVRI and ugriz ground-based light curves from the full monitoring campaign in AB magnitudes. Only the measurement uncertainties in the differential fluxes are shown. These uncertainties include those due to intercalibration, summarized in Table LABEL:tab:caldat. Systematic uncertainties for the absolute flux calibration are given in Table LABEL:tab:lcprop.
Table 2: WC18 ref ref ref ref LT ref ref ref ref ref LCOGT1 0.9 0.2 1.1 1.2 0.5 0.1 0.3 LCOGT2 2.3 0.4 1.8 0.5 1.9 0.4 1.9 LCOGT3 1.5 0.6 LCOGT4 2.7 0.3 0.1 0.2 LCOGT5 0.6 1.3 0.7 LCOGT6 0.5 1.6 1.4 0.4 LCOGT7 1.9 LCOGT8 1.3 WMO 0.9 0.6 0.3 CrAO 0.2 0.3 0.3 RCT 0.2 0.1 MO15 0.4 0.7 FWO 0.3 HLCO 0.4 KAIT 0.1 MLO 0.7 LOAO 0.4 \@ifx@empty\@ifx@empty Note. — Percentages are averaged for all epochs of the given telescope, measured relative to the flux at that epoch—see §2.3 for the definition of the intercalibration uncertainty. “ref” is the reference telescope to which the others are aligned. Note. — Percentages are averaged for all epochs of the given telescope, measured relative to the flux at that epoch—see §2.3 for the definition of the intercalibration uncertainty. “ref” is the reference telescope to which the others are aligned.
Table 3: u 56684.78 LT -34438.0 1292.5 56685.79 LT -24994.0 1050.5 56686.77 LT -29223.0 707.5 B 56645.64 WC18 -16159.0 460.55 56646.65 WC18 -16055.0 764.35 56647.65 WC18 -22095.0 595.81 g 56684.78 LT -8759.7 1384.1 56685.79 LT -8305.6 192.33 56686.77 LT -8892.6 579.47 V 56645.62 WC18 -4676.7 338.64 56646.61 WC18 -8876.5 500.65 56647.63 WC18 -6668.4 524.11 r 56684.78 LT -16832.0 491.06 56685.79 LT -23984.0 915.05 56686.77 LT -21698.0 1786.4 R 56644.64 CrAO -7835.7 4544.9 56646.63 WC18 -10378.0 588.18 56647.64 WC18 -13352.0 585.57 i 56684.78 LT -14639.0 2267.2 56685.78 LT -13055.0 2074.4 56686.77 LT -9307.6 458.99 I 56645.63 WC18 -6575.6 499.37 56646.63 WC18 -3043.8 732.31 56647.64 WC18 -5582.9 702.53 z 56684.78 LT -4697.2 984.0 56685.78 LT -3818.5 197.04 56686.77 LT -2659.7 1013.4 \@ifx@empty\@ifx@empty Note. — This table is available in its entirety in the online version of this article. Note. — This table is available in its entirety in the online version of this article.
Table 4: 1157.5 56690.61 56691.54 56692.39 1367. 56690.61 56691.54 56692.39 1478.5 56690.65 56691.58 56692.41 1746. 56690.65 56691.58 56692.41 \@ifx@empty\@ifx@empty Note. — This table is available in its entirety in the online version of this article. Note. — This table is available in its entirety in the online version of this article.

2.4 Light-Curve Properties

Table LABEL:tab:lcprop gives a summary of the sampling properties of the AGN STORM continuum light curves, and shows that the light curves have approximately daily cadence over the entire campaign. Paper I and Paper II only presented the HST 1367 Å continuum light curve; here, we include three additional UV continuum light curves measured from the HST COS spectra, extracted from 5–6 Å windows centered at 1158 Å, 1479 Å, and 1746 Å, and given in Table LABEL:tab:HST_lightcurves. These continuum windows were chosen to be as uncontaminated as possible by absorption lines and broad emission-line wings. We also drop the Swift V-band light curve from this analysis, because its mean fractional uncertainty is much larger than that of the ground-based Johnson V-band light curve (3.2% compared to 0.8%). The reported wavelengths of the optical light curves are pivot wavelengths calculated from the filter response curves of the optical bands, and they are independent of the source spectrum (atmospheric cutoffs at 3000 Å and 1 m were imposed for these calculations). Figure 2 shows a comparison of all the continuum light curves used in this study.

Table LABEL:tab:lcprop also gives the variability properties of the light curves. Column 8 gives the mean flux and root-mean square (rms) scatter of the light curves, corrected for Galactic extinction assuming a Cardelli, Clayton, & Mathis (1989) extinction law with and  mag (Schlegel et al. 1998; Schlafly & Finkbeiner 2011; Paper I). Columns 10, 11, and 12 give different estimates of their fractional variability. The fractional variability of a light curve is defined by


and the uncertainty in is


where is the value of the light curve at epoch , is the associated uncertainty, is the (unweighted) mean value of the light curve, and is the mean square of the measurement uncertainties (Rodríguez-Pascual et al., 1997; Vaughan et al., 2003). We also estimated the fractional variability using the JAVELIN amplitudes, , since this is an equivalent measure of under the DRW model. The values of and are often in good agreement, but with notable exceptions, as given in Table LABEL:tab:lcprop.

Figure 3 shows the mean flux and variability properties of these light curves. The top panel displays the mean SED (corrected for Galactic extinction). The vertical error bars show the minimum and maximum states of the AGN, which occur at and , respectively. These dates are based on the HST 1367 Å light curve, and the other bands are adjusted for interband time delays that are measured in §3. The middle panel illustrates the logarithm of the difference in flux between the minimum and maximum states of the AGN, which cleanly isolates the variable component of the spectrum and better traces the shape of the accretion-disk SED. For comparison, a standard thin accretion disk SED with is shown, arbitrarily normalized to match the Johnson V-band differential flux. Although the data are in excellent agreement with the prediction at longer wavelengths, the UV data lie significantly below the model SED. This discrepancy may be caused by extinction internal to the AGN, or the inner edge of the disk, which will display an exponential Wien cutoff rather than . A more complete discussion and modeling of the variable spectrum will be presented by Starkey et al. (in prep.). Finally, the bottom panel shows as a function of wavelength. The far-UV light curves have values of , which sharply decrease with wavelength to about 0.06 in the V band. At longer wavelengths, the trend flattens, reaching 0.02 in the z band.

At least part of this effect is caused by the constant flux contributed by the host galaxy, which becomes increasingly important at longer wavelengths. Based on spectral decomposition models and synthetic photometry (described in §4.1 and §4.2), the host galaxy contributes about 20% of the observed flux in the B band, and about 54% in the I and z bands. We corrected for this constant component, and Figure 3 and Table LABEL:tab:lcprop also show the host-galaxy flux and revised values of . The effect on the trend in Figure 3 is fairly subtle, and does not change the flattening at optical wavelengths.

The larger variability amplitudes at short wavelengths suggest that the SED of NGC 5548 becomes bluer in higher flux states. The same effect was seen by Cackett et al. (2015) in NGC 5548 with near-UV grism monitoring data from Swift. However, the trend is driven by the light curves at wavelengths  Å, and is most significant at wavelengths  Å, which may be why optical studies of AGN variability do not always find any “bluer when brighter” trend (e.g., Sakata et al., 2010).

Table 5: HST 1158 0.050 171 1.03 1.00 HST 1367 0.050 171 1.03 1.00 HST 1479 0.050 171 1.03 1.00 HST 1746 0.058 171 1.03 1.00 Swift UVW2 1928 0.030 284 0.62 0.39 Swift UVM2 2246 0.030 256 0.55 0.35 Swift UVW1 2600 0.030 270 0.52 0.38 Swift U 3467 0.020 145 1.20 0.99 Ground u 3472 0.035 270 0.52 0.38 Ground B 4369 0.030 151 1.11 0.98 Swift B Swift 4392 0.016 271 0.52 0.37 Ground g 4776 0.034 172 1.01 0.97 Ground V 5404 0.030 429 0.41 0.31 Ground r 6176 0.032 172 1.01 0.93 Ground R 6440 0.030 136 1.28 0.96 Ground i 7648 0.021 178 0.98 0.96 Ground I 8561 0.030 98 1.73 1.02 Ground z 9157 0.011 186 0.93 0.91 \@ifx@empty\@ifx@empty Note. — gives the number of epochs in the lightcurve, gives the average cadence, gives the median cadence, gives the mean flux (the uncertainty gives the rms scatter of the lightcurve), “Host” gives the host-galaxy flux, is defined in §2.4, and is the DRW amplitude. The flux calibration uncertainty is the systematic uncertainty for conversion to physical units (i.e., zeropoint errors). For HST, these values are taken from Paper I, while for Swift the values are from Table 6 of Poole et al. (2008). The uncertainties for the ground-based lightcurves represent our calibration to the SDSS AB mag photometric system. A correction for Galactic extinction has been applied to these data (see §2.4 for details). Note. — gives the number of epochs in the lightcurve, gives the average cadence, gives the median cadence, gives the mean flux (the uncertainty gives the rms scatter of the lightcurve), “Host” gives the host-galaxy flux, is defined in §2.4, and is the DRW amplitude. The flux calibration uncertainty is the systematic uncertainty for conversion to physical units (i.e., zeropoint errors). For HST, these values are taken from Paper I, while for Swift the values are from Table 6 of Poole et al. (2008). The uncertainties for the ground-based lightcurves represent our calibration to the SDSS AB mag photometric system. A correction for Galactic extinction has been applied to these data (see §2.4 for details). a Å b Corrected for host-galaxy flux
Figure 2. : AGN STORM UV and optical continuum light curves used in this analysis, restricted to the observing window of the HST campaign. Light curves have been converted to AB magnitudes, but are rescaled and shifted for clarity—the scales along the vertical axis show the fractional variations. The vertical dashed lines mark local extrema in the HST 1158 Å light curve.
Figure 3. : Top Panel: Mean SED of NGC 5548 from far-UV to optical wavelengths, corrected for Galactic extinction. The vertical error bars represent the AGN in the maximum and minimum states of the campaign. The horizontal error bars represent the rms width of the filter transmission curves. See §4 for a discussion of the host-galaxy estimate. Middle Panel: Variable SED component, calculated from the difference in flux between the minimum and maximum states, which more cleanly identifies the accretion-disk spectrum. The dashed red line is the predicted spectrum for a standard thin disk—discrepancies at short wavelengths may be due to extinction internal to the AGN or the inner edge of the disk. Bottom Panel: Fractional variability as a function of wavelength. is the DRW amplitude from the JAVELIN fits. For clarity, a small shift in wavelength to the points has been applied.

3 Time-Series Analysis

We measure the lags between light curves using two methods. First, we use the interpolated cross-correlation function (ICCF), as employed by Peterson et al. (2004), and estimate the uncertainty of the lag using a Monte Carlo method. Second, we use JAVELIN, which measures lags by modeling reverberating light curves as shifted, scaled, and smoothed versions of the driving light curve.

In the first method, the ICCF is calculated by shifting one light curve on a grid of lags spaced by 0.1 day, and calculating the correlation coefficient by linearly interpolating the second light curve. The lags are estimated from the centroid of the ICCF, defined as the mean ICCF-weighted lag for which . Uncertainties are estimated using the flux randomization/random subset selection (FR/RSS) method, wherein a distribution of ICCF centroids is built from cross correlating realizations of both light curves. Each realization consists of randomly selected epochs (chosen with replacement), and the corresponding flux measurements are adjusted by random Gaussian deviates scaled to the measurement uncertainties. The lags reported here correspond to the medians of the ICCF centroid distributions, while the lower and upper uncertainties define their central 68% confidence intervals.

We detrended the light curves, as is common practice in RM studies (Peterson et al. 2004; Paper II), in order to remove long-term secular trends that are poorly sampled in the frequency domain and may bias the observed lag (Welsh, 1999). The detrending procedure consists of subtracting a second-order polynomial linear least-squares fit (with equal weight given to all data points) from the observed light curves. Following Paper I and Paper II, we restricted the analysis to the time period coincident with the HST campaign, and measured the time delays relative to the HST 1367 Å light curve. When calculating the ICCF, we only interpolate the 1367 Å light curve. Table LABEL:tab:lags summarizes the resulting mean lags, corrected for cosmological time dilation (the redshift of NGC 5548 is ; Paper I). Lags for the hard and soft X-ray bands of the Swift XRT are also included, as determined in Paper II. The ICCF for all bands is shown in Figure 4 with the solid black lines, while the centroid distributions are shown as the gray histograms. We found that the HST 1158 Å and 1479 Å lags were only slightly larger than the spacing of our interpolation grid (0.1 day), so we repeated the procedure on these light curves using a grid of 0.01 day. This did not have a noticeable effect on the ICCF centroids, but it did change the ICCF peaks by day. The lags reported in Table LABEL:tab:lags make use of the finer grid for these light curves.

Our treatment of the Swift light curves (UVW2, UVM2, UVW1, U, and B) results in lags systematically larger than those found in Paper II, although the tension is only moderate (typically ). These differences are primarily caused by the different detrending procedures of the two studies. Paper II detrended the Swift light curves using a 30-day running mean, while we use a low-order polynomial. A running mean corresponds to a lower-pass filter than a polynomial trend, so this detrending procedure removes more low-frequency power from the light curve and is therefore expected to result in smaller lags. However, several of our light curves have very irregular sampling, which makes the calculation of the running mean poorly defined, so we instead use the low-order polynomial. The ground-based SDSS u and Swift U lags and the Johnson B and Swift B lags are consistent at the level using the polynomial detrending, so it is likely that the detrending procedure accounts for most of the difference between the near-UV lags. Two other smaller effects may be important for the lag determinations. First, the Swift UVOT optical filters are much narrower than the standard Bessell filters used for the ground-based light curves, so the observed variations are not perfectly identical (the Swift light curves also have slightly shorter baselines). Second, the Swift optical light curves have much larger fractional uncertainties, which may shift the ICCF centroid distribution of the otherwise similar light curves.

We also estimate the lags using JAVELIN, which calculates a maximum-likelihood lag, scale factor, and kernel width (assuming a top-hat transfer function) from the DRW covariance matrices. JAVELIN internally employs a linear detrending procedure, so we do not apply the second-order detrending as for the ICCF analysis. We also imposed a minimum kernel width of 0.75 day, in order to suppress solutions where JAVELIN finds a -function transfer function and aligns the reverberating light curve with the gaps between samples of the driving light curve (this is an aliasing problem associated with light curves that have similar cadences).

We adopt the medians of the posterior lag distributions and their central 68% confidence intervals as estimates of the lag and its uncertainty, which are given in Table LABEL:tab:lags. The posterior distributions are shown by the red histograms in Figure 4. The median lags are always consistent with the ICCF analysis, with the largest discrepancy being 1.7 in the r band. The Javelin uncertainties generically appear to be uncomfortably small. This is because JAVELIN assumes correctly characterized random Gaussian measurement errors, that the line light curve is a simple lagged and smoothed version of the continuum light curve, and that the smoothing kernel is well characterized by the functional form of the model (a top-hat function). Given that all these requirements are seldom fully met (particularly the Gaussianity of the measurement errors), Javelin uncertainties need to be interpreted conservatively. A rough rule of thumb from modeling gravitational lens time delays is that repeated measurements for the same system will typically be within 23 of each other.

The very small JAVELIN uncertainties may also indicate that the simple lagged and smoothed model of the reverberating light-curve model is an inadequate description of the data. Paper I found a similar result, where the shape of the line light curves was not always a good match to the observed continuum light curve. Therefore, smoothing the continuum light curve by a simple transfer function cannot always reproduce the line light curve, suggesting that other processes are important for the observed line emission (perhaps, for example, anisotropic emission/reprocessing). A more detailed investigation of this result will be pursued in upcoming papers of this series.

Using either lag estimation technique, we find that longer wavelength continuum variations follow those at shorter wavelengths. Figure 5 shows the lags as a function of the pivot wavelength of each filter. While the far-UV and near-UV light curves have time delays day, the V band lags the 1367 Å continuum by  days, and the z band lags it by days. For comparison, the He ii UV and optical lines (1640 Å and 4686 Å, respectively) have a mean lag of days relative to the 1367 Å light curve (Paper I, Paper V). The optical light curves have a time delay comparable to, and frequently larger than, that of high-ionization-state lines in the BLR.

The trend of larger lags at longer wavelengths is nearly monotonic. The most notable exceptions are (1) in the longest-wavelength filters, where the trend appears to level out near the i band, and (2) in the u and U bands. The u and U bands have mean lags of days and days, respectively, comparable to or larger than the lags of the g and V-band light curves. This may be due to emission originating in the BLR picked up in the u and U-band filters, which would contaminate measurements of the AGN continuum emission and artificially increase the observed lag. A similar explanation may exist for the downturn at the I and z bands, since Paschen continuum emission from the BLR begins at 8204 Å (see Korista & Goad 2001). We return to the question of BLR contamination in §4.

Optical continuum lags in NGC 5548 have previously been measured by Sergeev et al. (2005), and the same light curves were reexamined by Chelouche & Zucker (2013) and Chelouche (2013). Sergeev et al. (2005) found substantially longer time delays between the B and R/Cousins I bands than the lags presented here (about 8 days). However, the Sergeev et al. (2005) light curves have  day cadence and suffer from large seasonal/scheduling gaps of 20 days or more. The difference in the optical lags is therefore likely caused by systematic issues with the Sergeev et al. (2005) light curves, such as unfortunate gaps that affect the cross-correlation functions. On the other hand, Chelouche & Zucker (2013) and Chelouche (2013) claim that the large optical lags are due to BLR contamination and that the true continuum lags are consistent with zero. We discuss this possibility further in §4.3, but we find this interpretation to be unlikely. These studies did not discuss the impact of gaps in the data on the multivariate cross-correlation function used to disentangle line and continuum lags, and we are further skeptical that this method can meaningfully measure lags below the cadence of the light curves (3 days, in this case).

To avoid the systematics associated with small lags, interband continuum lags should be measured with data taken near or well below the timescale of any suspected lags. The UV wavelength coverage of the STORM project therefore lends a tremendous boost to our ability to detect the continuum lags, since the UV-optical lags are 3 to 6 times larger than the interband optical lags. This has implications for ground-based studies attempting to resolve interband continuum lags. Since the optical lags are of order 1 day (or less), the diurnal cycle may make it impossible to measure reliable interband optical lags from the ground without favorable conditions.

In order to quantify the trend of lag with wavelength, we fit a model to the data presented in Figure 5 using the functional form


where is the observed lag, is a reference wavelength, and and are free parameters. As in Paper II, we set  Å and report all covariances between parameters. The results of the fits are summarized in Table LABEL:tab:fitparams. We find that day and . The parameters are strongly correlated, with a normalized correlation coefficient , and , which approaches a low probability for 18 degrees of freedom ( and for a one-tailed test). Since there is good reason to suspect that the u and U bands are affected by BLR emission (see §4), we also fit the data excluding these lags. With these bands excluded, we find   day and  . The normalized correlation coefficient does not change, but the goodness-of-fit is now with , and (and for the same one-tailed test). The interpretation of Equation 4 is discussed in §5.3.

Figure 4. : ICCF for all light curves, with the ordinate showing the correlation coefficient . Lags for data from Paper I and Paper II (following our reanalysis) are shown in the left column, ground-based optical lags are presented in the right column. The grey histograms are the ICCF centroid distributions from the FR/RSS method, the red histograms are from JAVELIN. Both histograms are in units of .
Figure 5. : Time delay (ICCF centroid) as a function of pivot wavelength of the filters. The horizontal error bars represent the rms width of the filters. The best-fit model is shown by the dashed magenta line, while the fit fixing is shown by the dotted magenta line. Predictions for a thin-disk model with are shown by the solid cyan lines, although the assumptions of the model are unlikely to hold at large (see §5.3). The mean lag of the He ii and lines is shown by the horizontal dashed black line (Paper I, Paper V).
Table 6: Swift HX 4.4 Swift SX 25.3 HST 1158 HST 1479 HST 1746 Swift UVW2 1928 Swift UVM2 2246 Swift UVW1 2600 Swift U 3467 Ground u 3472 Ground B 4369 Swift B 4392 Ground g 4776 Ground V 5404 Ground r 6176 Ground R 6440 Ground i 7648 Ground I 8561 Ground z 9157 \@ifx@empty\@ifx@empty Note. — Measured relative to the HST 1367 Å light curve and corrected to the rest frame. The Swift lags are recalculated from Paper II using a second-order polynomial detrending routine, as described in §3. Note. — Measured relative to the HST 1367 Å light curve and corrected to the rest frame. The Swift lags are recalculated from Paper II using a second-order polynomial detrending routine, as described in §3.
Table 7: All -0.99 25.94 1.44 4/3 38.66 2.03 No Uu -0.99 16.85 1.05 4/3 22.64 1.33 No UuIR -0.99 12.4 0.89 4/3 13.0 0.87 \@ifx@empty\@ifx@empty

4 Contamination by Broad-Line Region Emission

As noted above, the u and U lags are outliers from the trend in Figure 5. A major component of the flux observed in these filters is the “small blue bump,” caused by bound-bound and bound-free hydrogen emission (the so-called Balmer continuum), as well as blended Fe ii lines that originate in the BLR. This BLR emission may cause the u and U-band lags to be biased estimators of the light-crossing time within the continuum source. In fact, several filters pick up other spectral features that originate in the BLR. The strongest is the prominent H line in the r and R bands, although additional emission lines and a diffuse continuum consisting of bound-free, free-free, electron scattering, and reflection is expected to be present at all wavelengths (Korista & Goad, 2001). Understanding the impact of BLR emission on the observed lags is therefore important for interpreting the interband time delays.

In this section, we assess the effect of BLR emission on the interband continuum lags. First, we decompose spectra of NGC 5548 into models of each emission component. We then estimate the fractional contribution from BLR emission in each filter using synthetic photometry. Finally, we simulate broad-band filter observations by combining mock continuum and BLR light curves, and search for biases in the lags by cross-correlating each emission component with the 1367 Å light curve.

4.1 Spectral Decomposition

We begin by decomposing spectra of NGC 5548 into models of each emission component. We obtained moderate-resolution () optical spectra of NGC 5548 using the Multi-Object Double Spectrographs (MODS; Pogge et al., 2010) on the Large Binocular Telescope (LBT; Hill et al. 2010). These observations are from 2014 June 08 and 2014 June 25 UT (HJD =2,456,817 and 2,456,834, respectively). The spectra were reduced and flux-calibrated using the modsIDL Spectral Reduction Pipeline.111A full description can be found at The spectra cover the wavelength range from 3000 Å to 1 m. Wavelength solutions were derived from comparison-lamp calibrations for each observing run. Relative-flux calibration was performed using three standard stars observed on the same nights as NGC 5548; however, the observations were taken in poor atmospheric conditions, making their absolute-flux calibration unreliable. We therefore rescaled the spectra so that the integrated [O iii]  fluxes match the value measured for the photometric nights of the optical spectral RM campaign, (Paper V). The slit width and extraction window of the MODS spectra were 5 and 15, respectively, chosen to match those of the optical monitoring spectra. This ensures that the relative contribution of host-galaxy light, narrow-line emission, and BLR emission are the same in both datasets. We corrected for Galactic extinction following the prescription described in §2.4. We did not make any correction for telluric absorption because broad-band filters suffer from the same effect.

Since we are only concerned with the relative magnitude of various emission components to the broad-band filter fluxes, we employed a minimal spectral decomposition, which is relatively coarse compared to state-of-the-art spectral modeling. Accordingly, we do not interpret any of our model parameters as indicative of physical conditions within the AGN, and instead focus on finding a model that provides a good fit to the data (based on minimizing ). Our decomposition has three components: host-galaxy starlight, the underlying AGN continuum, and the Balmer continuum shortward of  Å (rest frame). We ignore the diffuse continuum at other wavelengths, since it is poorly constrained, while the Balmer continuum can be determined from the shape and amplitude of the small blue bump. Emission-line fluxes are then estimated by subtracting the summed model components from the observed spectrum.

We simultaneously fit each component with an MCMC calculation, masking AGN emission lines and telluric absorption. We also masked the long and short edges of the spectra, because the MODS flux calibration is unreliable at  Å and  Å (rest frame). At these wavelengths, we set the observed flux equal to the summed model, which implicitly sets the emission-line flux to zero. This has a small effect on the estimated BLR contamination in the u, U, I, and z bands, but is more robust than using the unreliable flux calibration.

Details of the model components are as follows.

  1. Host Galaxy: We determined the host-galaxy spectrum using the STARLIGHT spectral synthesis code (Cid Fernandes et al., 2004). STARLIGHT fits the observed spectrum with a linear combination of a large library of synthetic stellar populations that span a wide range of ages and metallicities (150 templates from Bruzual & Charlot 2003). The best-fitting models consist of several very old (usually year) stellar populations at a range of metallicities (), and provide a reasonable match to the galaxy templates used by Denney et al. (2010) and Mehdipour et al. (2015). The resulting host templates have one parameter, the flux normalization. We also impose a tight prior on the flux at 5100 Å (rest frame), chosen to match the value measured by Bentz et al. (2013) adjusted to the MODS slit width and extraction window, (.

  2. Power Law: A broken power law is used to model the AGN continuum emission. This component has four free parameters: a flux-normalization factor, two spectral indices, and the location of the transition between indices. A loose prior (a Gaussian distribution with mean 5700 Å and width 700 Å) is imposed on the transition wavelength, to prevent it from moving to the edges of the spectra.

  3. Balmer Continuum: The Balmer continuum component is estimated from a grid of models calculated by Dietrich et al. (2002), evaluated at varying temperatures, electron densities, and optical depths. Again, we simply choose the template that produces the overall minimum value of . The templates have a single parameter, a flux rescaling factor.

We ignored blended Fe ii emission, because Fe emission is relatively weak in NGC 5548 (Denney et al., 2009; Mehdipour et al., 2015) and varies with an amplitude 50–75% that of H (Vestergaard & Peterson, 2005). This component is therefore expected to contribute very little flux to the broad-band photometric measurements and have a negligible impact on the observed lag. In order to assess the effect of this omission, we also fit the spectra with the small blue bump template of Mehdipour et al. (2015), which includes blended Fe ii emission lines. We found that these templates produce a poorer fit than the Dietrich et al. (2002) templates at the blue end of the spectrum, which may be a result of the limited wavelength coverage of our MODS spectra in the near-UV.

Each epoch was fit independently, and the resulting component parameters are in reasonable agreement, after allowing for the intrinsic variability of the power-law and Balmer continuum. The flux rescaling factors of the power-law continuum and galaxy templates are degenerate, so the prior imposed on the host-galaxy flux at 5100 Å (rest frame) does the most to constrain these parameters. Figure 6 shows an example of the decomposition, using the spectrum from 2014 June 08, overlaid with the filter transmission curves.

Figure 6. : Decompositions of the MODS spectra from 2014 June 08, showing the contribution of the model components to different filters. “BC” is the Balmer continuum, “PL” is the power law, “Host” is the host-galaxy component, and “Lines” are the AGN emission lines. The emission lines are estimated by subtracting the total model from the observed spectrum. Johnson/Cousins optical filter transmission curves (and Swift U) are shown by the dashed black lines, SDSS filters are shown by the dot-dashed lines. The Swift U and u bands are truncated at 3000 Å and the I and z bands are truncated at 1 m, in order to represent the atmospheric transmission cutoff.

4.2 Synthetic Photometry

Next, we estimate the contribution of each model component to the observed flux in each broad-band filter. We first reapply Galactic reddening to the model components, since differential extinction may affect the integrated flux across broad-band filters. We then calculate the observed flux using the synphot IRAF task and filter transmission curves for the calibration telescopes (WC18 BVRI filters and LT ugriz filters), truncated at 3000 Å and 1m to represent the atmospheric transmission cutoff. Uncertainties in the broad-band fluxes of individual components were estimated by resampling the posterior distributions of the model component parameters and rerunning synphot times.

Table LABEL:tab:fracs shows the results of our synthetic photometry. The “Total” column was calculated from the original spectrum, and the fractional contributions of individual components are reported relative to this value. The uncertainties represent the central 68% confidence interval of the resampled synthetic photometry distributions. The uncertainties are generally less than 1% because of the tight prior on the 5100 Å host-galaxy flux, which forces the galaxy template to be nearly constant and limits the variation of the other model components.

We do not consider effects of changing detector sensitivity with wavelength, since quantum efficiency curves for different instruments are usually much more variable than their filter transmission curves. Quantum efficiency will have the largest impact on the I and z filters, limiting the response of these filters at wavelength shorter than the cutoff imposed at 1 m. We investigated this effect by truncating the filter response at 9000 Å and repeating the experiment (essentially simulating a very steep quantum efficiency curve). We found that the final fractional contributions of the host/power-law components in these bands changes by 1% or less, and is therefore of minimal importance for our conclusions.

We find that the power-law component is dominant from the u band through the V band (% of the flux), although the host galaxy makes considerable contributions even in the B band (%). At longer wavelengths, the power-law component and host galaxy contribute roughly equal amounts of flux, except for the r and R bands, which include a substantial contribution from the H line: 20% in the r band and 15% in the R band. Line emission in all other filters is %. Balmer continuum emission accounts for about 19% of the flux in the u and U filters. The Mehdipour et al. (2015) blended Fe templates contribute % of the observed flux in the g, V, and r bands, confirming that Fe emission is a negligible component of the broad-band fluxes of this object.

4.3 Impact on Time Delays

The final step is to estimate the impact of BLR emission on the recovered interband time delays. First, we simulate light curves for the AGN continuum, Balmer continuum, and BLR emission models. We then sum the component light curves to reproduce light curves as would be observed in a given filter, and calculate the lag between the composite light curve and the HST 1367 Å continuum light curve.

The observed light curve is a superposition of the continuum emission and BLR emission,


where is the observed light curve in filter , is the continuum light curve in that filter, and is the line light curve, assumed to originate in the BLR. We use the term “line light curve” to refer to any emission produced in the BLR, including the Balmer continuum.

To simulate , we calculated the lag implied by the best-fit parameters in Figure 5 (  and   in Equation 4) at the pivot wavelength of the filter, and shifted the JAVELIN DRW model of the HST 1367 Å light curve by this amount. This method assumes that the HST 1367 Å light curve drives through instantaneous reprocessing after some light-travel-time delay, as would be expected for X-ray reprocessing in the accretion disk. 222 Reprocessed emission is also expected to be somewhat smoothed in time compared to the driving light curve. We therefore also considered versions of which are both smoothed and shifted by convolving the JAVELIN 1367 Å model with a top-hat function of amplitude for . We found that this smoothing made very little difference on the results, and so we only discuss the results for the shifted versions of here.

In RM, the line emission is assumed to be powered by ionizing continuum emission, so that


where is the driving continuum light curve and is the transfer function. For simplicity, we assume equal to the JAVELIN model of the HST 1367 Å light curve and a top-hat transfer function,


where is the mean line lag and is the width of the smoothing kernel. The choice of a top-hat function is for mathematical convenience and does not reflect any particular geometry, although it is widely consistent with a range of BLR configurations (for example, a spherical shell or the gross properties of an inclined disk/annulus; Peterson 2001). We varied and by octaves, with = 2, 4, 8, and 16 days, and = 0, 2, 4, 8, and 16 days. These values were chosen to sample the parameter space near the mean H lag during the monitoring campaign ( days; Paper V). To a low approximation, the Balmer continuum and H lag would be expected to lie near this value. Finally, we enforced causality by setting for .

We simulate light curves for the u, U, r, and R bands, in order to investigate the impact of Balmer continuum and H emission on the recovered lags. After generating the grid of shifted and smoothed line light curves, we renormalized each so as to reproduce the level of BLR contamination inferred from the spectral decomposition (Table LABEL:tab:fracs). We then adjusted the fractional variability amplitude (defined in §2.4) of both the continuum and line light curves to match their observed values. For the continuum light curves, is estimated directly from the observed broad-band light curves (Table LABEL:tab:lcprop, Column 9). For the line light curves, we set , derived from the observed H light curve (Paper V). We also experimented with changing the fractional variability amplitude of the line light curve to 0.012, 0.023, 0.092, and 0.184. Examples of two composite light curves and their model components, and , are shown in Figure 7.

After constructing and for each model, we calculated the lags of these light curves relative to using the ICCF method described in §3. In all cases, we recovered the input values of and to within the time resolution of the model light curves (0.12 day). We then calculated the ICCF for the composite light curve , finding that the recovered lags are most sensitive to the choice of and but are virtually independent of . The resulting mean lags are shown in Figure 8 as a function of input for the three values of near that of H (larger or smaller values of do not plausibly reproduce the observed lags, and are omitted for clarity). Larger values of these parameters tend to increase the recovered lag, but at the fiducial values of H the change is 0.6–1.2 days (blue point in Figure 8 with ).

We also checked for an effect of BLR contamination on the lag uncertainties. For each model, we found that larger values of and tend to increase the width of the ICCF. However, there was no correlation between these parameters and the location or width of the ICCF centroid distribution. This means that the lag uncertainties depend more sensitively on the light-curve quality rather than on any BLR contamination.

For values of that are smaller than ( in r and R and in u and U), the input line lag only has a limited effect on the recovered lag, evidenced by the flattening of the trends in Figure 8. This result is in contrast to the simple expectation that the observed lag is the flux-weighted mean lag of the line and continuum light curves, which scales linearly with the line lag. Instead, it appears that the observed lag only follows the line lag if the BLR emission dominates the variability properties of the composite light curve, as seen for the mock r and R bands at large . This indicates that the bias of the continuum lag introduced by BLR emission will usually be limited for broad-band filter light curves that are dominated by continuum emission, although the bias may still be important for small continuum lags.

Our simulations with these fiducial H parameters produce u and U-band lags in excellent agreement with the observed lags, while the simulated r and R-band lags overestimate the observed lag by about 1 day. Our current campaign cannot directly address the issue of the unknown values of and for these reverberations. However, it is expected from photo-ionization modeling that the Balmer continuum has a larger response () but shorter lag than H, while H should have a smaller response but longer lag (Korista & Goad, 2001, 2004). Based on Figure 8, this would serve to reduce the discrepancy between the recovered and observed lag in the r and R bands, while the recovered and observed lag in the u and U bands would remain in good agreement.

Thus, our simulations suggest that contamination by BLR emission can reasonably account for the systematic offset of the measured u, U, r, and R lags above the fit in Figure 5. This bias is well resolved in the u and U bands (the offset from the fit in Figure 5 is and , respectively), but of small importance in the r and R bands ( and , respectively). This result justifies our exclusion of the u and U-band data in the fit to Equation 4.

Chelouche & Zucker (2013) and Chelouche (2013) claim that BLR emission is responsible for the large B-R/Cousins I lags in the Sergeev et al. (2005) NGC 5548 light curves, and they find optical continuum lags consistent with 0 days. This is at odds with our results, since the BLR biases would have to be  days. These studies use a variation of the ICCF method (the multivariate CCF) to disentangle line and continuum lags from emission observed in a single filter. As we have already noted, gaps in the Sergeev et al. (2005) data make cross-correlation functions that rely on interpolation unreliable. Furthermore, this bias would imply that line emission contributes 30–50% of the flux in the R and Cousins I bands, which is implausibly high based on both our spectral decompositions and the composite Seyfert 1 spectrum of Chelouche (2013).

Table 8: 2014 June 08 U 8.42 u 8.43 B 7.23 g 7.40 V 7.39 r 8.91 R 8.44 i 7.45 I 5.77 z 4.60 2014 June 25 U 8.28 u 8.29 B 7.02 g 7.21 V 7.29 r 8.93 R 8.47 i 7.53 I 5.91 z 4.73 2014 June 08 (Blended Fe) U 8.49 u 8.49 B 7.23 g 7.40 V 7.38 r 8.90 R 8.43 i 7.44 I 5.77 z 4.59 2014 June 25 (Blended Fe) U 8.41 u 8.39 B 7.02 g 7.21 V 7.27 r 8.88 R 8.40 i 7.46 I 5.82 z 4.65 \@ifx@empty\@ifx@empty Note. — PL is power law, BC is Balmer continuum, Host is the host galaxy, Lines are AGN emission lines. BC includes a Fe emission template in the “Blended Fe” models. Note. — PL is power law, BC is Balmer continuum, Host is the host galaxy, Lines are AGN emission lines. BC includes a Fe emission template in the “Blended Fe” models.
Figure 7. : Examples of mock light curves, , , and , used for the analysis in §4.3. The top panel shows the HST 1367 Å light curve and the JAVELIN model used to generate the mock light curves, with the 1 uncertainty shown by the grey band. The middle panel displays an example of a mock u-band light curve, with a large line lag and high fractional variability, likely to result in the largest change of the observed lag. The bottom panel shows an example of a mock R-band light curve, with a more realistic line lag and fractional variability, chosen to be consistent with the H light curve. See §4.3 for further details.
Figure 8. : Recovered lags of mock light curves as a function of input line light curve lag. The colored points show the results for different variability amplitudes of the line light curve. The solid blue lines indicate the variability amplitude observed in the H light curve. The black dashed line represents the input lag of the continuum light curve , while the red dashed line is the observed lag and the red band is its 1 uncertainty. See §4.3 for further details.

5 Discussion

5.1 UV/Optical Light Curves and Lags

A primary goal of the AGN STORM project was to investigate how the continuum emission changes as a function of wavelength, and to assess any systematic issues introduced by using the optical continuum in place of the far-UV or extreme-UV. Figure 2 shows a detailed comparison of the HST 1367 Å light curve and all other data used in this study. We draw particular attention to the ground-based V-band light curve, since this is the most common choice of ionizing continuum proxy in ground-based RM studies. All of the major events and salient characteristics of the 1367 Å light curve are reproduced in the V band. There are, however, several noticeable differences.

5.1.1 UV–Optical lags

The first difference is a time delay between variations in the UV and optical light curves. Emission at 1158 Å, the shortest continuum wavelength available in this study, probably originates from a region of the accretion disk similar to that of the true ionizing continuum at  Å. This is because the lag-wavelength relation must flatten at small wavelengths (owing to the inner edge of the disk), but the inner edge already makes an important contribution to emission at  Å (Novikov & Thorne, 1973). Extrapolating the fit from Equation 4 to  Å implies a 0.26 day lag relative to the 1367 Å light curve, which is in reasonable agreement with the 1367 Å-1148 Å lag ( day). We therefore adopt a value of 0.2 day for the lag between the true ionizing continuum and the 1367 Å emission, since the lags for wavelengths  Å are unlikely to be much larger. This translates to a distance between the true ionizing continuum and the optically emitting portion of the disk of light days. A consequence of this UV-optical lag is that the radius of the BLR in NGC 5548 is underestimated when derived from the optical-H lag. The optical–H lag is variable in time, but typically has a value between 6 and 20 days (Peterson et al., 2004; Zu et al., 2011). Thus, if a similar UV–optical lag exists in other AGN, the physical size of the BLR is being systematically underestimated by up to (or for a lag of 20 days).

This result does not affect current optical RM SMBH masses, because RM only directly measures the virial product of the BLR, , where is the BLR lag and is its velocity dispersion (estimated from line-profile widths). Since the geometry and dynamics of the BLR are unknown, the virial product must be rescaled by a factor in order to produce a SMBH mass. While every AGN has a different value of , a statistical average can be calculated by calibrating an ensemble of virial products to some other SMBH mass estimate. Currently, this is done using the relation of local quiescent galaxies (Onken et al., 2004; Woo et al., 2010, 2013, 2015; Park et al., 2012a, b; Grier et al., 2013a). Thus, any systematic misestimation or bias of the lag (or velocity dispersion) is compensated by the calibration of , while the uncertainty of a single RM SMBH mass is dominated by the statistical uncertainty in , currently about (Grier et al., 2013a; Woo et al., 2015). However, any physical interpretation of (for example, a measure of the mean inclination of the BLR, assuming a disk or otherwise flattened geometry) requires a recalibrated value of that takes into account the UV–optical lag.

Single-epoch SMBH mass estimates are also unaffected by this result, since the radius-luminosity (RL) relation is inferred from a sample of RM AGN. While the larger BLR radius measured from the UV data would increase the normalization of the RL relation, a recalibration of exactly cancels this change. The UV–optical lag may introduce a second-order effect on single-epoch SMBH masses, if it is found that the magnitude of the UV–optical lag correlates with continuum luminosity or SMBH mass. Furthermore, the magnitude of the lag depends on accretion rate (see §5.3), which may also add scatter to existing mass-scaling relations. To investigate these effects, more simultaneous UV and optical RM experiments must be executed, using a sample of AGN with a wide range of luminosities.

Finally, the UV–optical lag has an impact on masses derived from direct dynamical modeling of RM data, since this method interprets the continuum-line lag as a measure of the physical radius of the BLR. To a low approximation, a larger BLR radius implies a proportionally larger SMBH mass. The effect of using UV continuum light curves for dynamical modeling studies will be investigated in future papers in this series, but until such modeling is complete, we adopt a RM-based SMBH mass for NGC 5548, since this estimate is less model dependent. From the H virial products compiled by Bentz & Katz (2015), and taking (Grier et al., 2013b), we adopt a mass of for the SMBH in NGC 5548. We note that this value moves in the correct direction for a larger BLR, but is still consistent within the quoted uncertainties of the dynamically modeled mass in Pancoast et al. (2014).

5.1.2 Optical Smoothing

The second difference between the UV and optical continuum light curves is that the V-band light curve appears to be smoother than the HST light curve. For example, the rapid oscillations in the UV light curve between HJD = 2,456,760 and 2,456,810 also appear in the V-band light curve, but at a much smaller amplitude with gentler inflections. The smoothing becomes increasingly severe at longer wavelengths where the amplitude of short-timescale variations decreases (see §2.4). These effects were also seen in NGC 2617 by Shappee et al. (2014), NGC 6814 by Troyer et al. (2016), and MCG-6-30-15 by Lira et al. (2015). Increased smoothing and decreased amplitudes are expected if shorter-wavelength emission drives the optical continuum, since the size, structure, and inclination of the accretion disk define a “continuum transfer function” that smooths the reprocessed light curve, while geometric dilution decreases the energy flux incident on large disk radii that contribute most to longer-wavelength emission.

In practical terms, the sharpest and strongest features in the V-band AGN STORM light curve are only slightly affected by this smoothing. Since these features provide the most leverage for constraining the CCF (Peterson, 1993), we conclude that the smoothing of the optical continuum is not important for ground-based RM studies that aim only to recover a mean emission-line lag and a SMBH mass. The smoothing may be more problematic for reconstructing velocity-delay maps, direct dynamical modeling, or regularized linear inversion (Horne et al., 1991, 2004; Bentz et al., 2010b; Grier et al., 2013b; Pancoast et al., 2014; Skielboe et al., 2015). These methods are very sensitive to the fine structure of the driving continuum light curve, and smoothing the light curve will erase information that would otherwise be helpful for reconstruction of the geometry and dynamics of the BLR. Velocity-delay maps, dynamical modeling, and regularized linear inversion for this dataset will be presented in upcoming papers in this series.

5.1.3 Magnitude of UV–Optical Lags

The large lags measured for optical bands, shown in Figure 5, are comparable to, and sometimes larger than, the lags for high-ionization-state lines such as He ii and C iv (Paper I). If the lags do in fact represent light-travel times across the accretion disk, then the optically emitting portion of the accretion disk appears to have a similar physical extent as the highly ionized portion of the BLR. This situation implies a close connection between the BLR and continuum-emitting source. For example, BLR clouds may be directly above or interior to the portion of the accretion disk emitting in the optical. Another plausible hypothesis is that at least part of the inner, high-ionization BLR emission arises from a wind launched from the surface of the accretion disk (e.g., Collin-Souffrin, 1987; Chiang & Murray, 1996; Proga & Kurosawa, 2010). Such models are able to reasonably explain many observed features of AGN emission lines, including their profiles, variability, and absorption characteristics (see Proga & Kallman, 2004; Eracleous et al., 2009; Denney, 2012; Higginbottom et al., 2014, and references therein). Alternatively, the accretion disk may smoothly merge with the BLR somewhere near 2–3 light days (for an analysis of this family of models, see, for example, Goad et al. 2012). Future papers in this series will attempt to map the geometry and kinematics of the inner BLR using the reverberation signal of high-ionization-state lines, which may shed further light on the connection between the accretion disk and BLR.

5.2 BLR Emission and Broad-Band Filter Lags

Based on our spectral decomposition, approximately 19% of the observed emission in the u and U bands is Balmer continuum emission from the BLR, while 15–20% of r and R-band emission is the prominent H line. These ratios may change with time, as shown in Table LABEL:tab:fracs, depending on the luminosity state of the AGN, the difference in phase between the continuum and line light curves, and the light curves’ variability amplitudes. For mean flux levels near the BLR contamination in the u, U, r, and R bands, as well as variability amplitudes and line lags that match the observed H light curve, our experiments with mock light curves indicate biases in the interband continuum lag of days.

These results depend on the assumption that all BLR emission light curves have properties similar to the H light curve. It is likely that the diffuse continuum actually has a stronger response but smaller lag than H, while H is expected to have a weaker response but larger lag (Korista & Goad, 2001, 2004; Bentz et al., 2010a). Since these parameters have offsetting effects, it is unlikely that the lag biases caused by BLR contamination are larger than the fiducial estimates presented here (see Figure 8). Future RM programs can test this result by specifically targeting the diffuse continuum and H emission, putting stronger constraints on their variability amplitudes and mean lags.

The systematic tendency for the u, U, r, and R band lags to sit above the fit in Figure 5 can therefore reasonably be explained by BLR contamination. In the case of the u and U bands, the offset from the fit to Equation 4 is large compared to the predicted lag (as well as the observational uncertainty), which supports our decision to exclude these data from the final model. On the other hand, the r and R-band offsets are much smaller, so the BLR bias probably makes little difference for our final model. Extending this reasoning to the B, g, and V-band filters, the BLR contamination is less than 10%, which would result in even smaller biases.

It is therefore unlikely that there are any important biases of the continuum lags in these bands, unless the diffuse continuum component (e.g., free-free emission or the Paschen continuum) makes a substantial contribution. This diffuse continuum component of the spectrum is unconstrained in our spectral decomposition, but it provides an intriguing possibility of explaining the downturn of the lag-wavelength relation in the I and z bands. The Paschen continuum begins at 8204 Å, between the i and I bands, so the true continuum lag-wavelength relation may run through the UV and Iz-band lags, but underneath the lags of the other optical filters. The viability of this explanation requires significant contamination of the optical filters by diffuse BLR emission, which can potentially be estimated through photoionization modeling of the HST data or additional optical/near-IR observations.

5.3 Accretion-Disk Size

A geometrically thin, optically thick, irradiated accretion disk makes definite predictions about the observed lag-wavelength structure of the AGN. Here, we compare this model to the observed continuum lags, although we do not necessarily interpret the model parameters as indicative of physical conditions within the AGN. Full physical modeling of the AGN STORM data is deferred to future papers in this series (Starkey et al., in prep.; Kochanek et al., in prep.).

The disk is assumed to have a fixed aspect ratio with scale height much smaller than radius, and is heated internally by viscous dissipation and externally by a UV/X-ray source near the SMBH at a small height above the disk. In such a scenario, the temperature profile is


where is the mass of the central SMBH, is the mass accretion rate of the disk, is the distance away from the black hole and central source of heating radiation, is the luminosity of the heating radiation, and is the albedo of the disk (Cackett et al., 2007). Here, we have ignored the inclination and the inner edge of the disk, as well as any relativistic effects. Inclination and relativity may have a small impact on the temperature profile, but the largest effect is caused by the inner edge, which reaches a maximum temperature and probably makes important contributions to emission at wavelengths  Å (Novikov & Thorne, 1973). This introduces an error when comparing the HST lags to this model, although the effect is small relative to the UV–optical lags.

Identifying the temperature with a characteristic emission wavelength , where is a multiplicative factor of order unity, and the radius with the light-travel time , we have


The factor accounts for systematic issues in the conversion of to for a given , since a range of radii contributes to emission at . From the flux-weighted mean radius


we derive , where is the inner edge of the disk, is the Planck function, and is the temperature profile defined in Equation 8.333Alternative definitions of exist. For example, a weighting function that better characterizes the radius responding to variable irradiation would replace Equation 8 with , and set equal to a constant fractional temperature variation. This yields .

If we measure relative to a reference time delay of a light curve with effective wavelength , then this becomes


Therefore, the parameter in Equation 4 is related to the energy generation rate responsible for heating the disk, while is predicted to be . The absolute size of the disk at can be measured by determining , which is inferred by assuming the corona is located at and fitting the X-rays lags (in which case ).

We can only determine indirectly through an estimate of the bolometric luminosity. We set , where is the radiative efficiency for converting rest mass into radiation, and quantifies all emergent radiation from the AGN, including coronal X-rays (in this sense, our model differs from the typical Shakura & Sunyaev 1973 thin-disk model). A convenient parameterization of is the Eddington ratio, . We also simplify Equation 5.3 by taking , where is the local ratio of external to internal heating, assumed to be constant with radius.

The equation for is then


A common choice for is 0.1, and we further assume that and for our fiducial calculations (i.e., the X-rays and viscous heating contribute equal amounts of energy to the disk). For a SMBH, these assumptions give day. If we increase the accretion rate by setting and 10, then and 0.65 day, respectively. The lag-wavelength relation for these models is shown in Figure 5, and the curves for –10 bracket our fit to Equation 4 with   days and  . However, it is important to note that the disk probably does not remain geometrically thin at these high accretion rates, and the assumptions of the model do not hold in this regime (Jiang et al., 2014; Sa̧dowski et al., 2014). Equation 12 is relatively insensitive to the ratio of external to internal heating—even if the X-rays contribute a negligible portion of the luminosity (), would only change by a factor of .

In Paper II, we found day, somewhat smaller than in this study. The smaller value can be explained by correlations between and , shown in Figure 9. For the final analysis in Paper II, was fixed to 4/3, and, if we do the same, we find day, in good agreement with Paper II. The fit with fixed has , making the lower value of  statistically preferred. However, this result is driven by the flattening of the lags at the reddest wavelengths. If we exclude the I and z bands from the fit (as well as u and U), we find and  day with , while fixing gives  day and .

Figure 9. : Probability distribution for the parameters and of our best-fit model. The contours show the 68, 95, and 99% confidence regions.

We therefore conclude that a reprocessing model can fit the data reasonably well but requires a much larger disk radius than predicted by standard thin-disk models. Fixing (in order to match the theoretical temperature profile), our best-fit value of day is a factor of 3.0 larger than the standard prediction with .

A sufficiently high accretion rate can account for this difference by increasing the size of the accretion disk. We note that uncertainties in the SMBH mass do not require to be larger than one, since (Equation 12), while the SMBH mass may be up to 1.75 times larger at 3 than our adopted value. This would still require to be somewhere in the range –1. On the other hand, a comparison of can be made assuming a thin-disk spectrum and using the observed optical luminosity (e.g., Collin et al., 2002; Netzer, 2013, Equation 4.53). From our spectral decompositions, we estimate that