Continuum Reverberation Mapping of Mrk 876 over Three Years with Remote Robotic Observatories

It is now known that most, if not all, galaxies each host a supermassive black hole (SMBH) in their respective nuclei (Richstone et al. 1998). An active galactic nucleus (AGN) occurs when material falls onto this SMBH. The conservation of angular momentum is expected to lead to the formation of an accretion disk. Most of these disks are too far away to be spatially resolved directly with current instruments, so indirect techniques are needed in order to observe and understand the AGN system. By measuring how the irradiated areas surrounding the SMBH respond to changes from the ionizing source, one can use time lags between different wavelength bands to measure the size scale. This technique is called reverberation mapping (Blandford & McKee 1982; Peterson 1993, 2014) and has been successfully carried out on a number of different AGN. Many regions of the AGN can be probed in this way, from the innermost areas of the accretion disk out to the dusty torus. For more information on reverberation mapping, a review can be found in Cackett et al. (2021).

Time lags between the continuum at different wavelengths are expected for reverberation of the accretion disk. In the lamppost model (Nayakshin & Kazanas 2002; Collin et al. 2003; Cackett et al. 2007), X-rays from an ionizing source located above the SMBH irradiate the accretion disk, where they are reprocessed and re-emitted at longer wavelengths depending on the temperature of the disk where they land. The disk is hotter closer to the SMBH, producing strong ultraviolet and continuum emission, while further and cooler regions of the disk predominantly produce the optical continuum. This model predicts that the X-rays should drive and lead the variability seen at longer wavelengths. Inner regions of the accretion disk should then see variations first, followed by regions further away from the SMBH. The measured time lag is then assumed to be dominated by the light travel time between different regions of the disk. A geometrically thin, optically thick accretion disk (Shakura & Sunyaev 1973) is expected to have a disk radial temperature profile of T(R) ∝ R^−3/4. Because wavelength maps to temperature via the Wien displacement law and time lag to radius, this gives wavelength-dependent lags, τ(λ), which should scale with wavelength as τ(λ) ∝ λ^4/3 (Collier et al. 1999). This relationship has been observed in previous continuum reverberation mapping studies, both in focused multi-instrument campaigns (e.g., Edelson et al. 2015; Fausnaugh et al. 2016; Hernández Santisteban et al. 2020; Kara et al. 2021) and from AGN survey studies (Jiang et al. 2017; Mudd et al. 2018; Guo et al. 2022). However, using a sample of over 9000 quasars from the SDSS Southern Survey Weaver & Horne (2022) find a bluer spectrum, τ(λ) ∝ λ^5/7, when possible Small Magellanic Cloud–like dust is considered.

Recent studies have shown some predictions of the lamppost model are not always valid, at least not for all AGN. Many AGN show a weaker correlation, or even a decoupling, between the X-ray and the ultraviolet/optical light curves (e.g., Edelson et al. 2017, 2019; Morales et al. 2019; Kara et al. 2021). Some AGN have trends indicative of systems more complicated than the lamppost model. Several studies have shown that there appears to be an incoming slow-moving lag, where the optical leads the ultraviolet and X-rays over the span of 100–1000 s of days (e.g., Breedt et al. 2009; Hernández Santisteban et al. 2020; Neustadt & Kochanek 2022). This implies that, on longer timescales, it is changes in the accretion disk's accretion flow that drive the variability, not an ionizing source. There is also evidence for obscuring elements, such as disk winds, that could exist between the ionizing source and the accretion disk, severing their direct connection (e.g., Kara et al. 2021). Previous reverberation mapping campaigns have found lags that indicate accretion disks ∼3 times larger than expected (e.g., Shappee et al. 2014; Edelson et al. 2015; Fausnaugh et al. 2016; Jiang et al. 2017). However, certain physical models have been created (Kammoun et al. 2019, 2021a, 2021b) that give lags consistent with observations. More studies involving long-term multiwavelength campaigns spanning several years are needed in order to understand variability on different timescales in AGN.

Mrk 876 is a prime target for such a study. Its location in the sky allows it to be observed by northern ground-based robotic observatories nearly year-round. It has a history of reliable variability and has been studied in several other independent campaigns (Kaspi et al. 2000; Bentz et al. 2013; Jha et al. 2022; Landt et al. 2023). From 2016 to 2019, it was observed by the Las Cumbres Observatory, the Liverpool Telescope, and the Dan Zowada Memorial Observatory. By combining these campaigns, a long-term look into the variability patterns and reverberation lags in Mrk 876 can be performed. In Section 2, we describe the data reduction and analysis. In Section 3, we present the results from time lag and spectral variability analysis. In Section 4, we discuss the implications of our findings. We summarize our results in Section 5.

Observations of Mrk 876 were taken from 2016 March through 2019 May. The observatories involved with this project are the Dan Zowada Memorial Observatory (Zowada), the Liverpool Telescope (LT), and the Las Cumbres Observatory (LCO). All images were processed using the standard pipelines for each observatory, which include bias/dark subtraction and flat-fielding. Each observatory uses the SDSS ugri filters. Zowada and LCO each have a Pan-STARRS z_s filter, while LT uses an SDSS z filter. Between two and five images per filter were taken each night, depending on the instrument, filter, and weather conditions.

Zowada is a robotic 0.5 m telescope located outside of Rodeo, New Mexico (Carr et al. 2022). It is owned and operated by Wayne State University. Located in the Canary Islands, LT is a robotic 2 m telescope operated by the Astrophysics Research Institute of Liverpool John Moores University (Steele et al. 2004). LCO is a global network of robotically operated telescopes (Brown et al. 2013). For this project, observations were taken from the 2 m telescope at Haleakala Observatory (ogg02), as well as the two 1 m telescopes (elp08 and elp06) at McDonald Observatory. Monitoring of Mrk 876 continued past 2019 May as part of a larger coordinated reverberation mapping campaign, and therefore those data are not used here for time lag analysis—they are only used to improve intercalibration between the light curves from the different telescopes. The LCO data from 2016 to 2018 were also used in Landt et al. (2023) to analyze the infrared reverberation signal. A summary of the observatories and observation epochs can be found in Table 1.

The light curves are obtained using differential photometry. A selection of comparison stars are chosen to be contrasted with the brightness of Mrk 876. We assume that the combined flux from the comparison stars is constant over time. We create an AGN light curve by calculating the flux of the AGN relative to the total flux of the comparison stars, allowing us to create the light curve for the AGN. Comparison stars are typically chosen such that they are 2–5 times brighter than the AGN in order to maximize the signal-to-noise ratio. All of the stars must be present in each telescope's field of view, meaning the choice of comparison stars was limited. Three stars were used for differential photometry based on the above factors. Different comparison stars are chosen for the u band, since most stars are typically fainter in this band.

For each image, the stars are found using the photutils (Bradley et al. 2021) module DAOStarFinder. Once found, we identify the comparison stars and AGN from the list of detected objects. A circular annulus and aperture are created for each object. The sizes of the apertures and annuli vary depending on the telescope. The annuli are chosen to have an inner radius of 20 pixels and an outer radius of 30 pixels for Zowada images. LCO and Liverpool have these annuli adjusted to match the same angular size based on their respective pixel scales. To choose the aperture size, we measure the signal-to-noise ratio for different aperture sizes and the fractional standard deviation of the comparison stars. We combine this in quadrature with the statistical uncertainty in flux. The aperture size with the lowest total flux uncertainty is chosen for each band and averaged to be used for each observatory. For Zowada, this is 5 pixels, for LT this is 8 pixels, for the LCO 1 m observatories (elp06 and elp08) this is 11 pixels, and for the LCO 2 m observatory (ogg02) this is 9 pixels. The median background is measured within the annulus and is scaled to the area within the aperture and subtracted. The average count rate is calculated from all of the observations of a specific filter taken on a given night, and these average observations are collected into a light curve for each band and each telescope.

We combine and intercalibrate the light curves from all the telescopes using CALI (Li et al. 2014). CALI assumes the AGN variability is described by a dampened random walk process in order to interpolate between gaps in data and align multiple telescopes' data to a common scale. It applies both additive and multiplicative factors to the data to achieve this. A Bayesian framework with a diffusive nested sampling algorithm is used to determine the intercalibration factors. Additional systematic errors may exist, so CALI increases the uncertainty on all measurements with a systematic error term that is added in quadrature to the original uncertainties. The complete set of light curves can be found in the left panel of Figure 1. When combined, we get an average cadence of 3.31 days in the g band. When we exclude large observational gaps of >30 days, we have an average cadence of 2.53 days.

Zoom In Zoom Out Reset image size

3.1. Measured Time Lags

The light curves for the ugriz filters can be found in Figure 1. To quantify the variability, the excess variance F_var (Edelson et al. 2002; Vaughan et al. 2003) is calculated for each of the light curves. Strong variability is observed, with respective variability amplitudes of 13%, 19%, 16%, 13%, and 12% in the u, g, r, i, and z bands. To determine the time lags, the cross-correlation function (CCF) is found for these light curves using the g band. This band is chosen as the comparison band because it has the best-sampled light curve and the highest variability amplitude. The CCF is found using the Python module PyCCF (Sun et al. 2018), which follows the Peterson et al. (2004) implementation of the ICCF technique to determine lag uncertainties. PyCCF calculates the CCF of two unevenly sampled light curves using interpolation to fill in the gaps between data points. The mean of the CCF at 80% of the CCF's peak (the centroid) is taken as the time lag between the two respective bands. We find in Figure 1 that the CCF of Mrk 876 is extremely broad, with a nearly flat top spanning several hundred days. This indicates the lag is dominated by long-term variations observed in the light curve. The broad CCF prevents a robust reverberation lag measurement. The data are detrended by subtracting a moving boxcar average to remove this long-term trend. The process of detrending AGN for improved short-term lag measurements is a common practice (e.g., Welsh 1999; Peterson et al. 2004; McHardy et al. 2014; Edelson et al. 2015; Fausnaugh et al. 2016; Hernández Santisteban et al. 2020; Pahari et al. 2020). The data were boxcar subtracted with a variety of boxcar widths ranging from 25 to 250 days, in intervals of 25 days. The lags of each set of detrended light curves were calculated using PyCCF as well. All boxcar widths between 75 and 175 days produced CCF lags consistent within 1σ uncertainties. The boxcar width with the lowest average lag uncertainty was 100 days, and therefore it was chosen for continued analysis. This process does not guarantee that the resulting time lags are true representations of the intrinsic lag between the different wavelengths. While the majority of boxcar lengths tested produce time lags that agree within 1σ, other detrending approaches may produce different results. The light curves and CCF of the 100 days detrended data can be found in Figure 2. A plot comparing all of the recovered time lags versus detrend lengths is shown in Figure 3, with the lags given in Table 2. The g lags are used for time-lag comparison against the other bands. We also measure the lag of the g band against itself, allowing the g band to have a lag beyond exactly 0 days. We use the detrended light curves for most analyses going forward. When we use the non-detrended data, we will refer to these as the base data. The light curves span several years, and due to Mrk 876's position in the sky there are annual seasonal gaps when Mrk 876 was not visible. The CCF method simply performs a linear interpolation between the gaps. However, more sophisticated methods, e.g., JAVELIN and PyROA, have been developed to use the variability properties of the light curves to inform a more realistic interpolation. JAVELIN (Zu et al. 2011, 2016) uses a dampened random walk model for the power spectrum of the light curves, and it assumes a top-hat transfer function. A recent comparison between CCF and JAVELIN methods for determining lags is presented by Yu et al. (2020). We use the Python 3 implementation of JAVELIN for our analysis. PyROA (Donnan et al. 2021) uses a running optimal average to determine a model for the driving light curve, and it assumes that each light curve is a time-shifted and flux-scaled version of this model. The Bayesian Information Criterion is then used to determine how much smoothing is required for the data, and model parameters are estimated using Markov Chain Monte Carlo. The time lags calculated using each method are given in Table 3, while in Figure 4 we show the lags as a function of wavelength. We also give the maximum correlation coefficient (${R}_{\max }$) in Table 3. All the light curves are well correlated. The lags generally agree between each method within their errors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Time Lag Comparison

Boxcar Width	Lag (days)
(Days)	u	g	r	i	z
0 (Base)	${14.60}_{-3.63}^{+3.72}$	${0.00}_{-2.52}^{+2.50}$	$-{5.17}_{-2.68}^{+2.85}$	${0.43}_{-2.52}^{+2.87}$	${2.31}_{-2.97}^{+3.10}$
25	$-{0.29}_{-0.80}^{+0.81}$	${0.00}_{-0.36}^{+0.35}$	${0.55}_{-0.48}^{+0.49}$	${2.18}_{-1.12}^{+0.92}$	${1.27}_{-1.31}^{+2.69}$
50	$-{1.79}_{-0.85}^{+0.73}$	${0.00}_{-0.44}^{+0.43}$	${2.07}_{-0.51}^{+0.48}$	${6.48}_{-0.90}^{+0.98}$	${6.51}_{-2.11}^{+1.49}$
75	$-{4.05}_{-0.83}^{+0.79}$	${0.00}_{-0.38}^{+0.39}$	${3.14}_{-0.45}^{+0.47}$	${10.03}_{-0.84}^{+0.89}$	${7.42}_{-1.27}^{+1.31}$
100	$-{5.52}_{-0.89}^{+0.87}$	${0.00}_{-0.34}^{+0.34}$	${3.32}_{-0.43}^{+0.42}$	${10.95}_{-0.88}^{+1.14}$	${7.68}_{-1.06}^{+1.20}$
125	$-{5.34}_{-1.01}^{+0.90}$	${0.00}_{-0.38}^{+0.39}$	${3.24}_{-0.44}^{+0.48}$	${11.6}_{-1.06}^{+1.33}$	${8.12}_{-1.12}^{+1.19}$
150	$-{4.11}_{-1.08}^{+1.05}$	${0.00}_{-0.41}^{+0.41}$	${3.26}_{-0.48}^{+0.49}$	${12.21}_{-1.05}^{+1.41}$	${9.24}_{-1.18}^{+1.51}$
175	$-{3.56}_{-1.16}^{+1.05}$	${0.00}_{-0.41}^{+0.40}$	${2.86}_{-0.47}^{+0.48}$	${11.03}_{-1.09}^{+1.34}$	${10.12}_{-1.55}^{+1.62}$
200	$-{2.26}_{-1.20}^{+1.28}$	${0.00}_{-0.46}^{+0.45}$	${2.04}_{-0.59}^{+0.62}$	${8.21}_{-1.07}^{+1.08}$	${8.87}_{-1.63}^{+1.43}$
225	$-{1.50}_{-1.48}^{+1.40}$	${0.00}_{-0.62}^{+0.63}$	${1.07}_{-0.83}^{+0.80}$	${6.74}_{-1.34}^{+1.15}$	${5.8}_{-1.63}^{+1.53}$
250	$-{0.43}_{-2.27}^{+2.20}$	${0.00}_{-1.01}^{+1.01}$	${0.02}_{-1.35}^{+1.25}$	${4.90}_{-1.94}^{+1.6}$	${2.80}_{-2.14}^{+1.81}$

Notes. A comparison of the PyCCF time lags for different detrending lengths. The first row with a detrend length of 0 days represents non-detrended data, which we will refer to as the base data in further analyses.

Download table as:  ASCIITypeset image

Table 3. Detrended Time Lags and Light-curve Properties

Method	u	g	r	i	z
${{\rm{R}}}_{\max }$	0.75	1.0	0.95	0.80	0.68
PyCCF	$-{5.52}_{-0.89}^{+0.87}$	$-{0.0}_{-0.34}^{+0.34}$	${3.32}_{-0.43}^{+0.42}$	${10.95}_{-0.88}^{+1.14}$	${7.68}_{-1.06}^{+1.20}$
JAVELIN	$-{1.66}_{-0.06}^{+3.96}$	$-{0.0}_{-0.00}^{+0.00}$	${3.51}_{-0.03}^{+0.99}$	${10.49}_{-2.03}^{+2.01}$	${7.5}_{-1.03}^{+3.97}$
PyROA	$-{4.86}_{-0.78}^{+0.78}$	$-{0.01}_{-0.16}^{+0.16}$	${3.18}_{-0.40}^{+0.25}$	${9.46}_{-0.69}^{+0.62}$	${8.63}_{-0.56}^{+1.15}$

Notes. R ${}_{\max }$ is the maximum value of the cross-correlation coefficient between the two light curves (calculated with respect to g). A value of 1 is maximal correlation, while a value of 0 means no correlation between the signals. The lags for PyCCF, JAVELIN, and PyROA are all given in units of days and are calculated with respect to the g band.

Download table as:  ASCIITypeset image

3.2. Theoretical Time Lags

The expected time lag τ can be calculated from a simple model of the accretion disk. Assuming a standard lamppost-like X-ray corona ionizing the accretion disk, the time lag should relate to the radius of the disk R as τ ∼ R/c. For a geometrically thin, optically thick accretion disk (Shakura & Sunyaev 1973), the temperature profile follows $T{(R)\propto (M\dot{M})}^{1/4}{R}^{-3/4}$. Assuming blackbody radiation such that λ ∝ T⁻¹, one finds a relationship between the time lag and the emitted wavelength of light:

$$\begin{eqnarray}&&\tau \propto {(M\dot{M})}^{1/3}{T}^{-4/3}\propto {(M\dot{M})}^{1/3}{\lambda }^{4/3}.\end{eqnarray} \tag{ 1 }$$

We test this relation by fitting the following function:

$$\begin{eqnarray}&&\tau ={\tau }_{0}[{(\lambda /{\lambda }_{0})}^{\beta }-{y}_{0}],\end{eqnarray} \tag{ 2 }$$

where τ₀ is the normalization parameter and is measured in days, β is the relationship of wavelength to the measured time lags, λ is the wavelength of the band being observed, and λ₀ is the reference wavelength band used, which for this study is λ₀ = 4770 Å (the effective wavelength of the g band). The adjustment factor y₀ is present to prevent the lags from being exactly zero at λ₀, with its value normally being around unity. A standard thin disk predicts β = 4/3. We fit the lag-wavelength relation both fixing β = 4/3 and allowing it to be a free parameter in the fit. The best fits are shown in Figure 4 and given in Table 4. When β is left as a free parameter, it drops below unity, leading to an unreasonably large value for τ₀. As such, it is not included in the figure. This equation holds the assumption that the X-ray emitting component's height is small relative to the radius of the accretion disk. If this is not true, then this acts as a lower limit to the measured lags.

Table 4. Fitted Accretion Disk Properties

Method	τ0 (days)	H (RG )	${\dot{m}}_{{\rm{E}}}$
PyCCF	9.22 ± 2.64	54.0 ± 56.0	0.416
JAVELIN	6.50 ± 2.00	28.0 ± 55.0	0.416
PyROA	9.16 ± 1.38	48.0 ± 25.0	0.416

Notes. Comparison of the fitted parameters between PyCCF, JAVELIN, and PyROA. The normalization parameter τ₀ from Equation (1) is found when β is fixed to be 4/3, as expected from thin disk theory. These are plotted as the solid black lines in Figure 4. The X-ray corona height (H) is found when the mass accretion fraction ($\dot{m}$) is fixed to be 0.416 for a spin parameter a* = 0.998. For more details, see Section 3.2.

Download table as:  ASCIITypeset image

We also fit to the lags the analytic prescription described by Kammoun et al. (2021b). This prescription predicts time lags using five different parameters. These are the black hole mass, mass accretion rate ($\dot{m}$), X-ray luminosity in the 2–10 keV range, X-ray corona height (H), and black hole spin (a*). We fit our observed time lags for when the spin parameter a* = 0.998 and a* = 0. The black hole mass is estimated to be 2.2 × 10⁸ M_⊙ (Peterson et al. 2004; Bentz & Katz 2015), and the X-ray luminosity in the 2–10 keV range is 10^44.11 erg s⁻¹, based on archival XMM-Newton data (Laha et al. 2018). We estimate the mass accretion rate using our calculation of the Eddington accretion rate (${\dot{m}}_{{\rm{E}}}$), which is described in Section 3.3 as a part of the analysis on the normalization parameter τ₀. The initial fits are done fixing $\dot{m}$ to be this value while H remains the only free parameter. This fit fails for a* = 0, with the corona height H becoming nonphysical. This agrees with previous investigations indicating that Mrk 876's SMBH has a high spin parameter (Bottacini 2022). We also perform the fitting allowing ${\dot{m}}_{{\rm{E}}}$ to be a free parameter. For this regime, we used a grid search to find the lowest possible χ² value. We search from 0 to 250 R_G for H and between 0 and 0.75 Eddington accretion for ${\dot{m}}_{{\rm{E}}}$. The upper and lower bounds for the uncertainty of the measurement are determined by finding the value of the parameters when χ² is 2.3 above the lowest value found. When left as free parameters, uncertainties on ${\dot{m}}_{{\rm{E}}}$ and H are largely unconstrained, but they agree within 1σ (for both zero spin and maximally spinning cases). We therefore only give the parameters when ${\dot{m}}_{{\rm{E}}}$ is fixed and a* = 0.998 in Table 4.

3.3. Calculation of Normalization Parameter τ₀

One noted problem among AGN reverberation mapping campaigns is the measured value of normalization parameter τ₀. We can estimate the expected value of τ₀ for the thin disk model using estimates of the black hole mass and mass accretion, and compare it to the measured value (Fausnaugh et al. 2016). The majority of campaigns have found that the fitted value of τ₀ is 2–3 times larger than the expected/calculated value (e.g., Shappee et al. 2014; Edelson et al. 2015; Fausnaugh et al. 2016; Jiang et al. 2017). This implies that the accretion disk itself is 2–3 times larger than the standard thin disk model predicts. Alternatively, some other aspect of the AGN system may be interfering with measurements and creating an erroneously large measurement of τ₀.

Fausnaugh et al. (2016) parameterize the following equation for calculating τ₀:

$$\begin{eqnarray}&&{\tau }_{0}=\displaystyle \frac{1}{c}{\left(X\displaystyle \frac{k{\lambda }_{0}}{{hc}}\right)}^{4/3}{\left[\left(\displaystyle \frac{{GM}}{8\pi \sigma }\right)\left(\displaystyle \frac{{L}_{\mathrm{Edd}}}{\eta {c}^{2}}\right)\left(3+\kappa \right){\dot{m}}_{{\rm{E}}}\right]}^{1/3}.\end{eqnarray} \tag{ 3 }$$

To calculate this value for Mrk 876, several assumptions are made. It is assumed that the X-rays and viscous heating contribute roughly the same amount of energy to the disk, such that the radiative efficiency for converting rest mass into radiation η = 0.1 and the local ratio of external heating to internal heating κ = 1. The factor X is a multiplicative factor of order unity, and it is determined from temperature T and the wavelength measured via T = Xhc/k λ. This factor helps account for how temperature relates to the wavelength emitted for a given radius, and it is influenced by the choice of geometry in the disk. For a flux-weighted mean radius, X = 2.49 (Fausnaugh et al. 2016), but taking variation of the disk emission into account leads to X = 5.04 (Tie & Kochanek 2018).

To calculate the Eddington accretion rate (${\dot{m}}_{{\rm{E}}}$), one must estimate the bolometric luminosity. We use the standard bolometric correction via L_bol ∼ 9λ L_λ (5100 Å) (Kaspi et al. 2000). However, because 5100 Å is not a central wavelength of the SDSS filters, the g band is used as the nearest available approximation. The g-band data are converted from relative flux to magnitudes using the comparison stars, which have magnitudes from the APASS catalog. The AGN magnitudes are then converted into fluxes and extinction corrected with an E(B − V) of 0.027 (Petrosian et al. 2007) using Cardelli's extinction law (Cardelli et al. 1989), and then corrected to rest-frame fluxes using z = 0.1385 (Lavaux & Hudson 2011). We assume a luminosity distance of D_L = 588.4 Mpc to calculate L(5100 Å) from the dereddened, rest-frame flux. The Eddington luminosity is calculated assuming a black hole mass of 2.18 × 10⁸ M_⊙ (Bentz & Katz 2015). For our bolometric luminosity, L_bol = 1.144 × 10⁴⁶ erg s⁻¹, we determine an Eddington fraction of ${\dot{m}}_{{\rm{E}}}=0.416$. Substituting into Equation (3), we calculate τ₀ = 2.57 days and 6.58 days for X = 2.49 and X = 5.04, respectively.

3.4. Spectral Analysis

We note that, in the time lag analysis the i-band lags (Figure 4) are consistently offset from the rest of the trend. At the redshift of Mrk 876, Hα is close to the effective wavelength of the i band. Other emission lines may be affecting the other lags as well. Figure 5 shows spectra taken on 2019 July 02 (just after the end of the campaign) from the LCO Haleakala Observatory (FTN) overlaid with the filters used to take the data. The z band shows no significant emission line contribution. The u, g, and r bands show some emission lines. To determine total flux contribution, we modeled the emission lines using astropy Spectrum1D models. These are shown in Figure 5 as the orange lines. The total contribution is summed up from each model for each point, using the throughput of each filter as a modifier of the total strength of the emission. The percentages of flux that come from the continuum versus the emission lines are found in Table 5. It should be noted that we only factor the broad emission lines into the emission line percentage, but not any potential contribution from the diffuse continuum, which is also thought to originate from the BLR. We find that the majority of filters see a small amount of emission line contribution, but not enough to warrant additional consideration. The exception to this is the i band, where we find that Hα contributes 33% of the total flux. To ensure that the presence of Hα was consistent throughout the campaign, we also analyze spectra taken from the start of the campaign in 2016. We find the Hα line contributes 29% of the total flux in 2016, confirming that Hα is a strong, consistent presence throughout the entire monitoring campaign.

Zoom In Zoom Out Reset image size

Table 5. Emission Line Contributions by Filter

Filter	Continuum	Emission Line
u	91%	9%
g	94%	6%
r	85%	15%
i	67%	33%

Note. The percentage of continuum and emission line contribution to the overall image from the modeled LCO spectra found in Figure 5.

Download table as:  ASCIITypeset image

However, a 30% contribution by Hα does not imply a 30% increase in lag. It is a common expectation that, to zeroth order, the continuum lag and the Hα lag will combine weighted by their flux; however, simulations have shown that F_var is the dominating factor (Fausnaugh et al. 2016). The Hα lag for Mrk 876 has been measured to be ${43}_{-22}^{+40}$ days, with a measured variability amplitude smaller than the continuum variability (Kaspi et al. 2000). However, it is more complicated than this, as the lag also depends on the variable flux, properties of the driving light curve, and the shapes of the transfer functions. Detailed simulations would be needed in order to properly assess this, and they are beyond the scope of this paper. Given that the Hα flux is a significant fraction of the flux in the i band, it is plausible to attribute the excess i-band lag to the Hα line.

3.5. Flux–Flux Analysis

To determine the spectral energy distribution (SED) of the variable flux, we perform a flux–flux analysis on the base data, similarly to Cackett et al. (2020) (see also Starkey et al. 2017; Hernández Santisteban et al. 2020). The photometric light curves are first flux-calibrated using the magnitudes of the comparison stars as found in the APASS catalog (Henden et al. 2018) for the g, r, and i bands, and the Pan-STARRs catalog (Chambers & Pan-STARRS Team 2018) for the z band. Neither catalog contained our u-band comparison stars, so as a proxy we use an observation from the Neil Gehrels Swift Observatory (Gehrels et al. 2004) of Mrk 876. The flux-calibrated light curves are corrected for Galactic absorption with an E(B − V) of 0.027 (Petrosian et al. 2007). We use the extinction law of Cardelli et al. (1989), and adjust the data to rest-frame flux. We then perform the flux–flux analysis by breaking the flux into constant and variable components representing the galaxy and the AGN, respectively, using the following formula:

$$\begin{eqnarray}&&{f}_{\nu }(\lambda ,t)={A}_{\nu }(\lambda )+{R}_{\nu }(\lambda )X(t).\end{eqnarray} \tag{ 4 }$$

A_ν is the average spectrum, R_ν is the rms spectrum, and X(t) is a dimensionless light curve normalized to a mean of zero and a standard deviation of unity. The light curves and fits for the non-detrended data are shown in panel (a) of Figure 6, and the flux–flux relations are shown in panel (b). To estimate the galaxy contribution to the different bands, we extrapolate the fits to where the uncertainty envelope of the shortest wavelength band crosses f_ν = 0, which we define as X(t) = X_G . This serves as a reference point for the other bands, and determining f_ν at X(t) = X_G provides a lower limit on the constant component in each band. The dashed lines X_F and X_B represent the lowest and highest points found from all filters, and X₀ is given as reference. Panel (c) of Figure 6 shows the maximum, minimum, and average SED of Mrk 876 during the monitoring, along with the variable (rms) and constant spectral components determined from the flux–flux analysis. Table 6 gives the values determined with the flux–flux analysis. The rms spectrum is consistent with the f_ν ∝ λ^−1/3 expected for an accretion disk spectrum. An excess in the variable spectrum in the i band is seen, which would be consistent with a significant variable broad Hα line contributing in that band. There is also an additional constant component in the i band shown in panel (c), indicating there is another source of continuum emission beyond the accretion disk.

Zoom In Zoom Out Reset image size

Table 6. Flux–Flux Results

	u	g	r	i	z
Max	5.58 ± 0.02	4.53 ± 0.01	4.79 ± 0.01	7.69 ± 0.02	5.74 ± 0.04
Mean	4.61 ± 0.00	3.52 ± 0.00	3.87 ± 0.00	6.56 ± 0.0	4.85 ± 0.00
Min	3.54 ± 0.02	2.22 ± 0.03	2.68 ± 0.03	4.9 ± 0.07	3.54 ± 0.09
Constant	0.09 ± 0.07	0.21 ± 0.01	0.78 ± 0.01	2.64 ± 0.02	2.08 ± 0.04
rms	0.92 ± 0.02	0.67 ± 0.00	0.63 ± 0.00	0.80 ± 0.01	0.56 ± 0.01

Notes. Base light curve flux–flux analysis values, which are shown in Figure 6, panel (c). All units are in milliJanskies (mJy).

Download table as:  ASCIITypeset image

We performed photometric monitoring of the AGN Mrk 876 over a 3 year period, during which it exhibited large-amplitude (12%–19%) variability in each of the ugriz bands (see light curves in Figure 1). We looked for lags between the different bands, as would be expected from an accretion disk reverberation scenario where ionizing radiation drives variability at longer wavelengths, with the hottest, inner (ultraviolet) part of the disk responding before the cooler, outer (optical) region of the disk.

Initial CCF lags recovered from PyCCF reveal a long-term variability (>100 ∼ days) that dominates over the short-term variability expected from accretion disk reverberation. AGN are known to vary on longer timescales, and to accurately recover the reverberation lags expected on shorter timescales (days), this long-term variability needs to be removed. We remove the long-term variability by subtracting a moving boxcar average from the light curves. To do this, a moving boxcar average is subtracted from the data. We test a range of widths to see which provided the lowest uncertainties in the resulting time lags. We find that detrend widths of 75–175 days produce time lags that agree within 1σ uncertainties, and we select 100 days because this width provides the lowest uncertainties. The base data CCF distributions are shown to the right of the light curves in Figure 1, and the detrended light curves and their CCFs can be found in Figure 2. Once the long-term trends are removed, the resulting CCF is significantly narrower, allowing a precise recovery of the reverberation lags. The measured lags from all tested detrending lengths can be found in Table 2.

The source of the long-term variability may be changes in the accretion rate. This acts on the viscous timescale, predicted to be on the scale of hundreds of days for the optical emitting region (Krolik et al. 1991). As more reverberation mapping studies are undertaken, it is becoming apparent that the lamppost model does not adequately explain all variability seen in some AGN light curves. It is possible that many AGN exhibit these long-term trends, as many observation campaigns detect some kind of long-term variation, often associated with accretion disk flow or broad-line region interference (Arévalo et al. 2008; Breedt et al. 2009; McHardy et al. 2014; Hernández Santisteban et al. 2020). It is difficult to perform the lengthy monitoring needed to capture these long-term variations with traditional observing on a single telescope. With the rise of robotic observatories like LT, LCO, and Zowada, as well as all-sky surveys such as ASAS-SN, ATLAS, ZTF, PAN-STAARS, and CRTS, more studies like this one will be possible.

The wavelength dependence of the detrended lags are shown in Figure 4. As a comparison, the lags are also found with the more sophisticated JAVELIN and PyROA techniques. We fit the relation τ ∝ λ^β to all sets of lags. The lags recovered by PyCCF, JAVELIN, and PyROA are generally in agreement within uncertainty, and all recovered lags are well represented by τ ∝ λ^4/3. They are consistent with the expected wavelength dependence for a Shakura–Sunyaev geometrically thin optically thick accretion disk with an illuminating central source. All methods find an excess in their i-band lags, deviating from an extrapolation of the trend in the other wave bands. The lags from different methods are given in Table 3.

Notably, we do not detect a u-band excess lag in our detrended time lag measurements. These excesses have been detected in a number of reverberation mapping studies (Edelson et al. 2015; Fausnaugh et al. 2016; Edelson et al. 2017; Cackett et al. 2021); however, they also have not been detected in some as well (Kara et al. 2021). Without UV/X-ray monitoring, we lack the wavelength coverage to determine if Mrk 876 truly lacks a u-band lag excess. However, if the source of the u-band excess is from the BLR (Korista & Goad 2001), then it is possible that our detrending process has already removed this contribution. Looking at the base light-curve time lags in Table 2, we see that u-band lag has the largest lag of all the bands at 14.68 days. This could imply that the BLR u-band emission operates on timescales of around 100 days, and that its contribution is removed by our detrending process. However, the large flat-topped CCFs that are produced prevent a robust lag measurement, so this result should be taken with caution.

The boundary of the accretion disk and the dusty torus can be examined using our data. We can extrapolate what the values for longer wavelength lags would be, assuming they follow the τ = λ^4/3 trend predicted for a geometrically thin optically thick disk. We compare the values we predict against the values determined from fitting the near-infrared rms spectrum by Landt et al. (2023). The implied radii are measured to be ∼25 light-days in the J band and ∼56 light-days in the K band. We extend our fits of τ ∝ λ^4/3 to the optical data to estimate the expected near-IR disk lags for these bands. The uncertainty on τ₀ and y₀ are used to create the upper and lower bounds for the lags. We extrapolate lags from each of the lag methods. For PyCCF, J = 17 ± 7 days and K = 48 ± 16 days. For JAVELIN, J = 14 ± 7 days and K = 36 ± 14 days. For PyROA, J = 17 ± 4 days and K = 47 ± 8 days. Our extrapolation of the disk lags is in general agreement with what is measured by Landt et al. (2023).

We also fit the time lags to the analytical prescription described in Kammoun et al. (2021b) in Figure 4. The black hole mass and X-ray luminosity can be found in the literature, leaving the mass accretion fraction (${\dot{m}}_{{\rm{E}}}$) and X-ray corona height (H) to be fit. We calculated the Eddington fraction using a bolometric luminosity ${\dot{m}}_{{\rm{E}}}$ = 0.416, so we fit for both when $\dot{m}$ is fixed to this value and when it is allowed to be a free parameter. However, when leaving ${\dot{m}}_{{\rm{E}}}$ as a free parameter, the fit is poorly constrained. We therefore only consider fits with ${\dot{m}}_{{\rm{E}}}$ fixed at 0.416. Figure 4 only shows the fit for when ${\dot{m}}_{{\rm{E}}}$ is a fixed parameter, with the red dashed line representing a* = 0.988 regime. All of the values determined from fitting can be found in Table 4. The values of H are found to be around 30–50 R_G , but they agree with each other to within 1σ.

Spectra of Mrk 876 (Figure 5) show that there is a significant contribution from the Hα broad line that is present in i-band measurements. This emission remains significant throughout the duration of the campaign. While determining the exact impact is beyond the scope of this paper, emission from the broad-line region has been suspected to influence lags in prior continuum reverberation mapping campaigns. Contributions from broad emission lines can also influence the lag in a photometric band (Chelouche 2013). While this is not prominent in all objects, the redshift of Mrk 876 puts the strong Hα line in the middle of the i band, indicating that it is a strong possibility in this case.

We perform a flux–flux analysis on the dereddened and redshift-corrected flux in Figure 6. The values determined for the SED are recorded in Table 6. The flux–flux analysis allows a determination of the variable and constant components of the SED. The variable (rms) component agrees with a f_ν ∝ λ^−1/3 spectrum expected for an accretion disk, though it shows an excess in the i band. This indicates the presence of additional variability beyond what is expected from the accretion disk. The constant spectrum also shows an excess in the i band. Both the variable and constant component excesses can be attributed to a prominent Hα line.

This excess variability lends credence to the broad Hα line being the source of the longer than expected i-band lag. The Hα line has variations smaller than what our detrending removes but longer than what the accretion disk is expected to create. Our detrending process removes slow variations on timescales around 100 days, while the Hα in Mrk 876 has been measured to vary on timescales of roughly ${43}_{-22}^{+40}$ days (Kaspi et al. 2000). This allows its variations to influence the i-band lags still, despite the detrending process. As explained in Section 3.4, we do not expect the lags of Hα and the accretion disk to add together simply, but the effect on the lags is clear to see.

We estimate the Eddington fraction using the bolometric luminosity to be ${\dot{m}}_{{\rm{E}}}=0.416$ during the campaign. This makes it one of the highest Eddington rate AGN studied via continuum reverberation mapping to date. Based on its black hole mass and mass accretion rate, like many other studies have found (Fausnaugh et al. 2016; Cackett et al. 2018; Edelson et al. 2019; Cackett et al. 2020; Kara et al. 2021), we find that the normalization (τ₀) recovered from fitting the measured time lags (9.22 ± 2.64 days, 6.50 ± 2.00 days, and 9.16 ± 1.38 days for PyCCF, JAVELIN, and PyROA, respectively) is several times larger than τ₀ = 2.57 days calculated from theory when X = 2.54. These values range from 2.6 to 3.6 times greater than theory, depending on the method. This implies that a flux-weighted mean radius alone cannot adequately describe the measured accretion disk sizes with the other assumptions about accretion disks we applied. In order for a flux-weighted mean radius model of the accretion disk to recover the size of the disk measured, we would need to assume a much higher accretion rate than expected, a higher Eddington ratio, or a lower accretion efficiency. This problem exists beyond just reverberation lag measurements, as similar results are found through gravitational microlensing campaigns (e.g., Poindexter et al. 2008; Mosquera et al. 2013).

The additional phenomena required are broadly divided into two categories. First are theories that the source of the lags is still X-ray reprocessing, but that some aspect of AGN geometry is either incorrectly assumed or different than expected. One example of this would be that the irradiating source is higher above the accretion disk than usually assumed (Kammoun et al. 2019). The other category of explanations involve different sources for the lags, such as lags instead being due to disk turbulence (Cai et al. 2020) or due to far-UV illumination from the inner disk shining onto and providing the reprocessing radiation for the rest of the disk (Gardner & Done 2017). Alternatively, Gaskell & Harrington (2018) suggest these longer-than-expected lags can be attributed to an underestimation of the intrinsic flux of the AGN, and hence an underestimation of the Eddington ratio, due to the presence of large amounts of intrinsic reddening. However, the fact that the variable spectrum closely follows f_ν ∝ λ^−1/3 would seem to suggest that there is not a large amount of intrinsic reddening in this particular object.

For our analysis in particular, the choice of detrending via a moving boxcar average could influence the lags and therefore the measured size of the accretion disk. While the majority of detrending lengths we test agree with each other within uncertainty (Figure 3), the outlier cases show that the different lengths produce smaller time lags. We argue that the best detrend length is that of 100 days, due to it having the smallest uncertainties, but again this does not guarantee that it is the correct length. The equation we use to calculate the lags (Equation (3)) is also simplistic in its assumption about the geometry of the AGN and accretion disk system.

Another possible explanation for the accretion disk size problem is an underestimation of X. Many studies use the value calculated by Fausnaugh et al. (2016) of 2.54. However, other studies have calculated it considering other factors of the AGN. When including variation of the disk emission (Tie & Kochanek 2018), the value of X then becomes 5.04. Using this value, we calculate τ₀ = 6.58 days. This is closer to what we observe for all lag measurement methods, falling within 1σ for all methods except PyROA.

Given its mass and Eddington fraction, the continuum lags in Mrk 876 are some of the longest yet observed. The outer edge of the accretion disk is expected to become self-gravitating at 12 light-days, regardless of the mass of the system (Lobban & King 2022). Our u to z lag is around 13 days, and the measured τ₀ of around 9 days suggests that the z band corresponds to a disk size of 17–20 days (depending on lag method), significantly larger than the 12 light days self-gravitating radius. Our estimates are consistent with what Landt et al. (2023) estimate through spectral fitting.

In summary, Mrk 876 displays large-amplitude variability over 3 yr and shows significant time lags across the optical band. The lags in Mrk 876 are some of the longest continuum lags yet measured, and longer than expected for the self-gravitating radius. Our conclusions on Mrk 876 are as follows:

• 1.  
  We measure broad CCFs with the base light curves, indicating long-term variations are present. The data are detrended by subtracting a moving boxcar average of 100 days, recovering the short-term lags. We measure the lags with PyCCF, JAVELIN, and PyROA. The results can be found in Table 3 and are plotted in Figure 4.
• 2.  
  The i-band lag is longer than expected from an extrapolation of the other bands. We analyze spectra taken from before and after the campaign, finding that the Hα broad-line emission has a strong contribution in the i band—up about 1/3 of the total i-band flux. We find that, due to its intermediately long lags (∼40 days), this signal would not be removed by detrending and would still exist in our lag measurements of the detrended data.
• 3.  
  We perform a flux–flux analysis on both the base and detrended data. The base data contains an i-band excess in both the constant and variable (rms) emissions, while the rest of the bands agree with the expected profile for an accreting thin disk f_ν ∝ λ^−1/3. The detrended flux–flux analysis reveals that the excess variable emission is removed via detrending, implying that what remains as an excess in the detrended light curves varies on scales longer than reverberation mapping but shorter than 100 days. This further supports the notion that Hα is responsible for this lag excess.
• 4.  
  The normalization parameter τ₀ is found for all lag measurement methods. We calculate this value following the parameterization described by Fausnaugh et al. (2016). Two different values of the factor X, 2.49 (Fausnaugh et al. 2016) and 5.04 (Tie & Kochanek 2018), are used to calculate τ₀ = 2.57 days and τ₀ = 6.58 days, respectively. Our τ₀ values are closer to the latter, agreeing for most methods within uncertainty.
• 5.  
  The lags are fit to the analytical prescription described in Kammoun et al. (2021b). When ${\dot{m}}_{e}$ is fixed at the observed value of 0.416, we find an X-ray source height of 30–50 R_G for a maximally spinning black hole.

Continuum reverberation mapping continues to challenge the standard picture of AGN accretion. More studies with high-cadence observations are required in order to truly understand the AGN system.

We thank the anonymous referee for providing comments and suggestions. J.A.M. and E.M.C. gratefully acknowledge support from the National Science Foundation through AST-1909199. We thank David Moutard for feedback that has improved this manuscript. This research made use of Photutils, an Astropy (Astropy Collaboration et al. 2013, 2018) package for detection and photometry of astronomical sources (Bradley et al. 2021). This work makes use of observations from the Las Cumbres Observatory global telescope network. Research at UC Irvine is supported by NSF grant AST-1907290. H.L. acknowledges a Daphne Jackson Fellowship sponsored by the Science and Technology Facilities Council (STFC), UK. E.R.C. acknowledges support by the NRF of South Africa. T.T. acknowledges support from the NSF through grant NSF-AST-1907208. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

Facilities: LCOGT (elp06 - , elp08 - , FTN - Faulkes Telescope North, optical) (Brown et al. 2013) - , Liverpool:2m (Steele et al. 2004) - Liverpool 2 meter telesope at Observatorio del Roque de los Muchachos, Swift (Gehrels et al. 2004) - Swift Gamma-Ray Burst Mission, Zowada (Carr et al. 2022) - .

Software: Astropy (Astropy Collaboration et al. 2013, 2018, 2022), CALI (Li et al. 2014), JAVELIN (Zu et al. 2011, 2016), photutils166 (Bradley et al. 2021), PyCCF (Sun et al. 2018), PyROA (Donnan et al. 2021).